Title: SAFEGUARDED LEARNED CONVEX OPTIMIZATION WITH GUARANTEED CONVERGENCE

Abstract: Many applications require quickly and repeatedly solving a certain type of optimization problem, each time with new (but similar) data. However, state-of-theart general-purpose optimization methods may converge too slowly for real-time use. This shortcoming is addressed by learning to optimize (L2O) schemes, which construct neural networks from parameterized forms of the update operations of general-purpose methods. Inferences by each network form solution estimates, and networks are trained to optimize these estimates for a particular distribution of data. This results in task-specific algorithms (e.g., LISTA, ALISTA, and D-LADMM) that can converge order(s) of magnitude faster than general-purpose counterparts. We provide the first general L2O convergence theory by wrapping all L2O schemes for convex optimization within a single framework. Existing L2O schemes form special cases, and we give a practical guide for applying our L2O framework to other problems. Using safeguarding, our theory proves, as the number of network layers increases, the distance between inferences and the solution set goes to zero, i.e., each cluster point is a solution. Our numerical examples demonstrate the efficacy of our approach for both existing and new L2O methods.

1. INTRODUCTION:

Solving scientific computing problems often requires application of efficient and scalable optimization algorithms. Despite the ever improving rates of convergence of state-of-the-art general purpose algorithms, their ability to apply to real-time applications is still limited due to the relatively large number of iterations that must be computed. To circumvent this shortcoming, a growing number of researchers use machine learning to develop task-specific algorithms from general-purpose algorithms. For example, inspired by the iterative shrinkage thresholding algorithm (ISTA) for solving the LASSO problem, a sparse coding problem, Gregor & LeCun (2010) proposed to learn the weights in the matrices of the ISTA updates that worked best for a given data set, rather than leave these parameters fixed. They then truncated the method to K iterations, making their Learned ISTA (LISTA) algorithm form a K-layer feed-forward neural network. Empirically, their examples showed roughly a 20-fold reduction in computational cost compared to the traditional algorithms. Several related works followed, also demonstrating numerical success (discussed below). While classic optimization results often provide worst-case convergence rates, limited theory exists pertaining to such instances of data drawn from a common distribution (e.g., data supported on a lowdimensional manifold). As a step toward providing such theory, this work addresses the question: Does there exist a universal method that encompasses all L2O algorithms and generates iterates that approach the solution set with guarantees?

We provide an affirmative answer to this question by prescribing and proving properties of neural networks generated within our L2O framework. Convergence is established by including any choice among several practical safeguarding procedures, including nonmonotone options. Nonmonotone safeguarding enables sequences to traverse portions of the underlying space where the objective function value may increase for a few successive iterations as long as, on average, the sequence approaches the solution set. Although counterintuitive, this ability may lead to faster convergence. Furthermore, we provide a practical guide in our discussion for how practitioners may use our framework to create and apply L2O schemes to their own problems.

The theoretical portion of this work is presented in the context of fixed point theory. This is done to be sufficiently general and provide the desired convergence result for the wide class of optimization methods that can be expressed as special cases of the Krasnosel'skiȋ-Mann (KM) method. For concreteness and ease of application, we then provide the special-case results to several well-known methods (e.g., proximal-gradient, Douglas-Rachford splitting, and ADMM).

Related Works.:

Learning to learn methods date back decades (e.g., see (Thrun & Pratt, 1998) for a survey of earlier works and references). A seminal L2O work in the context of sparse coding was by Gregor & LeCun (2010) . Numerous follow-up papers also demonstrated empirical success at constructing rapid regressors approximating iterative sparse solvers, compression, 0 encoding, combining sparse coding with clustering models, nonnegative matrix factorization, compressive sensing MRI, and other applications (Sprechmann et al., 2015; Wang et al., 2016a; b; c; d; Hershey et al., 2014; Yang et al., 2016) . A nice summary of unfolded optimization procedures for sparse recovery is given by Ablin et al. (2019) in Table A .1. However, the majority of L2O works pertain to sparse coding and provide limited theoretical results. Some works have interpreted LISTA in various ways to provide proofs of different convergence properties (Giryes et al., 2018; Moreau & Bruna, 2017) . Others have investigated structures related to LISTA (Xin et al., 2016; Blumensath & Davies, 2009; Borgerding et al., 2017; Metzler et al., 2017) , providing varying results dependent upon the assumptions made. Chen et al. (2018) introduced necessary conditions for the LISTA weight structure to asymptotically achieve a linear convergence rate. This was followed by Liu et al. (2019a) , which proved linear convergence of their ALISTA method for the LASSO problem and provided a result stating that, with high probability, the convergence rate of LISTA is at most linear. The mentioned results are useful, yet can require intricate assumptions and proofs specific to the relevant sparse coding problems.

L2O works have also taken other approaches. For example, the paper by Li & Malik (2016) used reinforcement learning with an objective function f and a stochastic policy π * that encodes the updates, which takes existing optimization algorithms as special cases. Our work is related to theirs (cf. Method 1 below and Algorithm 1 in that paper), with the distinction that we include safeguarding and work in the fixed point setting. The idea of Andrychowicz et al. (2016) is to use long short term memory (LSTM) units in recurrent neural networks (RNNs). Additional learning approaches have been applied in the discrete setting (Dai et al., 2018; Li et al., 2018; Bengio et al., 2018) . Balcan et al. (2019) reveal how many samples are needed for the average algorithm performance on the training set to generalize over the entire distribution. This is practical for choosing training data and may be used for training L2O networks within our framework.

Our Contribution. This is the first work to merge ideas from machine learning, safeguarded optimization, and fixed point theory into a general framework for incorporating data-driven updates into iterative convex optimization algorithms. In particular, given a collection of data and an update operator from an established method (e.g., ADMM or proximal gradient) for solving an optimization problem, we present procedures for creating a neural network that can be used to quickly infer solution estimates. The first novelty of this framework is the ability to incorporate several safeguarding procedures in a general setting. The second is that we present a procedure for utilizing machine learning methods to incorporate knowledge from particular data sets. However, our most significant contribution to the L2O literature is to combine these results into a single, general framework for use by practitioners on any convex optimization problem.

Outline. We first provide a brief overview of the fixed point setting of this work in Section 2. Then we present the SKM method and convergence results in Section 3. The incorporation of the SKM method into a neural network and subsequent training approach is presented in Section 4. This is followed in Section 5 by numerical examples, discussion in Section 6, and conclusions in Section 7.

2. FIXED POINT METHODS:

Let H be a finite dimensional Hilbert space (e.g., the Euclidean space R n ) with inner product •, • and norm • . Denote the set of fixed points of each operator T : H → H by Fix(T ) := {x ∈ H : T x = x}. In this work, for an operator T with a nonempty fixed point set (i.e., Fix(T ) = ∅), the primary problem considered is the fixed point problem: Find x ∈ Fix(T ).

(1)

Convex minimization problems, both constrained and unconstrained, may be equivalently rewritten as the problem (1) for an appropriate mapping T . The method chosen for solving the minimization problem determines the operator T in (1) (e.g., see Table 1 below for examples). We focus on the fixed point formulation to provide a general approach, given T , for creating a sequence that converges to a solution of (1) and, thus, also of the corresponding optimization problem.

The following definitions will be used in the sequel. A mapping T :

EQUATION

An operator T : H → H is α-averaged if α ∈ (0, 1) and there is a nonexpansive operator Q :

H → H such that T = (1 -α)Id + αQ, where Id is the identity operator. If the constant α is not important, then T may for brevity be called averaged. We also denote the distance between a point

x ∈ H and a set C by d C (x) := inf{ x -y : y ∈ C}.

(3)

Two operators frequently used in optimization are constructed from monotone relations. Letting α > 0 and f : H → R be a function, the resolvent of the (possibly) multi-valued subgradient ∂f is

J α∂f := (Id + α∂f ) -1 (4)

and the reflected resolvent of ∂f is

EQUATION

)

If f is closed, convex, and proper, then the resolvent is precisely the proximal operator, i.e.,

EQUATION

)

From these definitions, it can be shown that R α∂f is nonexpansive and J α∂f is averaged. Results above may all be found in Bauschke & Combettes (2017) (e.g., see Prop. 4.4, Thm. 20.25, Example 23.3, and Prop. 23.8) . See Table 1 for examples of these operators in optimization methods.

A classic theorem states that sequences generated by successively applying an averaged operator converges to a fixed point. This method comes from Krasnosel'skii (1955) and Mann (1953) , which yielded adoption of the name Krasnosel'skiȋ-Mann (KM) method. This result is stated below and can be found with various forms and proofs in many works (e.g., see (Bauschke & Combettes, 2017 , Thm. 5.14), (Byrne, 2008, Thm. 5 .2), (Cegielski, 2012, Thm. 3.5.4) , and (Reich, 1979, Thm. 2)). Theorem 2.1. If an averaged operator T : H → H has a nonempty fixed point set and a sequence {x k } k∈N with arbitrary x 1 ∈ H satisfies the update relation

EQUATION

then there is x ∈ Fix(T ) such that {x k } k∈N converges to x , i.e., x k → x .

There are pathological cases where the result fails for operators that are only nonexpansive (e.g., when x 1 = 0 and either T = -Id or T is a rotation). However, this is easily remedied since any convex combination of a nonexpansive operator with the identity is averaged.

3. SAFEGUARDED KM METHOD:

This section generalizes the classic KM iteration in (7). We accomplish this by defining an envelope of operators T L2O 

EQUATION

For ζ ∈ (0, 1), the operator T L2O (•; ζ) defined in ( 8) is ζ-averaged. Although using ζ 1 may result in an operator that fails to be averaged, this can be useful in accelerating the convergence of a method (e.g., see (Giselsson et al., 2016) ).

Table 1 : Below are the averaged operators for well-known algorithms. We assume α > 0 and, when α is multiplied by a gradient, we also assume α < 2/L, with L the Lipschitz constant for the gradient. The dual of a function is denoted by a superscript * and Ω := {(x, z) :

Ax+Bz = b}. The block matrix M is M = (α -1 Id, A T ; -A, β -1 Id). In each case, L is the associated Lagrangian. Problem Method Averaged Operator T min f (x) Gradient Descent Id -α∇f min f (x) Proximal Point prox αf min{g(x) : x ∈ C} Projected Gradient proj C • (Id -α∇g) min f (x) + g(x) Proximal Gradient prox αf • (Id -α∇g) min f (x) + g(x) Peaceman-Rachford R α∂f • R α∂g min f (x) + g(x) Douglas-Rachford 1 2 (Id + R α∂f • R α∂g ) min Ω f (x) + g(z) ADMM 1 2 Id + R αA∂f * (A T •) • R α(B∂g * (B T •)-b) minf (x) s.t. Ax = b Uzawa Id + α A∇f * (-A T •) -b minf (x) s.t. Ax = b Proximal Method J α∂L of Multipliers minf (x) + g(Ax) PDHG J M -1 ∂L

Method 1 Safeguarded Krasnosel'skiȋ-Mann (SKM)

1: Choose x 1 ∈ H and δ ∈ (0, 1)

2: for k = 1, 2, . . . do 3: Choose parameter ζ k 4: Choose µ k ∈ (0, ∞) 5: y k = T L2O (x k ; ζ k ) 6: if S(y k ) ≤ (1 -δ)µ k then 7:

x k+1 = y k 8:

else 9:

x k+1 = T (x k ) 10:

end if

11: end for Example 3.2. Let f : H → R be closed, convex and proper, and define

EQUATION

For fixed ζ ∈ (0, ∞), the operator

T L2O (•; ζ) is averaged.

The practicality of T L2O is discussed and demonstrated in Sections 4 and 5, respectively. In the remainder of this work, each operator T : H → H is assumed to be averaged and we set S := Id -T . Our proposed method below is called the Safeguarded Krasnosel'skiȋ-Mann (SKM) Method. Explanation of the SKM Method is as follows. In Line 1, the initial iterate and parameter δ are initialized. A common choice for the initial iterate is x 1 = 0. The for loop on Lines 2 to 11

Table 2 : Choices for µ k updates that ensure Assumption 3 holds. Here α ∈ (0, 1) and, for fixed m ∈ N, Ξ k is the set of the most recent min{m, k} indices for which the inequality in Line 6 holds.

NAME UPDATE FORMULA:

Geometric Sequence

µ k+1 = (1 -α)µ k if Line 6 holds, µ k otherwise. GS(α)

Decrease µ k by geometric factor whenever Line 6 holds.

Arithmetic Average

m k+1 := m k + 1 if Line 6 holds, m k otherwise.

AA

µ k+1 :=    1 m k + 1 S(x k+1 ) + m k µ k if Line 6 holds, µ k otherwise.

Use m k to count how many times Line 6 holds and µ k is the average of the residuals among those times.

Exponential Moving Average:

µ k+1 := α S(x k+1 ) + (1 -α)µ k-1 if Line 6 holds, µ k otherwise. EMA(α)

Average µ k with the latest residual whenever Line 6 holds.

Recent Term:

RT

µ k+1 = S(x k ) if Line 6 holds, µ k

otherwise. Take µ k to be most recent residual for which Line 6 holds.

Recent Max:

µ k+1 = max ∈Ξ k S(x )

RM(m):

Take µ k to be max of the most recent residuals for which Line 6 holds.

generates each update x k+1 in the sequence {x k } k∈N . The choice of parameter ζ k in Line 3 may be any value that results in a well-defined operator T L2O (•; ζ k ) in Line 5. The choice of µ k in Line 4 defines the safeguarding procedure that is used to ensure convergence. Safeguarding is implemented through a descent condition inequality in Line 6. When the inequality in Line 6 holds, T L2O (x k ; ζ k ) is used to update x k via Line 7. Otherwise, a KM update is used to update x k via Line 9. Notice also the iteration indexed parameters may all be chosen dynamically (rather than precomputed).

Below are several standard assumptions used to prove our convergence result in Theorem 3.1. Assumption 1. The operator T is nonexpansive with a nonempty fixed point set.

The following assumption ensures boundedness of sequences generated by the SKM method.

Assumption 2. The operator S is coercive, i.e.,

EQUATION

)

Remark 3.1. Assumption 2 does not hold, in general, for nonexpansive operators. For example, if T is the gradient operator (Id -α∇f ) for some α > 0 and f is a constant function, then S(x) = 0 for all x ∈ H. However, a minor perturbation to the f enables Assumption 2 to hold. In this example, if one fixes small ε > 0 and sets f (x) := f (x) + ε 2 x 2 , then the associated S satisfies S(x) = ε x . This idea generalizes and, since this works for arbitrarily small ε, in practice it may be reasonable to assume Assumption 2 holds when applying the SKM Method. Assumption 3. If the inequality in Line 6 is satisfied infinitely many times, then the sequence {µ k } k∈N converges to zero. Remark 3.2. Assumption 3 may be enforced by using various choices that are dependent upon combinations of the previous residuals. This is illustrated by Table 2 and Corollary 3.1 below.

Our main convergence result is Theorem 3.1 below (proven in the Appendix).

Theorem 3.1. If {x k } k∈N is a sequence generated by the SKM method and Assumptions 1 to 3 hold, then lim

EQUATION

And, if {x k } k∈N contains a single cluster point, then {x k } k∈N converges to a point x ∈ Fix(T ).

We propose several methods for choosing the sequence {µ k } k∈N in Table 2 . These methods are adaptive in the sense that each update to an iterate µ k depends upon the current iterate x k and (possibly) previous iterates. These update schemes enable each µ k to trail the value of the residual norm S(x k ) . This implies there may exist j ∈ N for which

EQUATION

)

Such leniency is desirable since it is possible that, even though the residual norms converge to zero, constructing a sequence of iterates along the quickest route to the solution set requires traversing portions of the underlying space where the residual norms increase for a few successive iterations.

The safeguarding schemes in Table 2 are justified by the Corollary below (proven in the Appendix).

Corollary 3.1. If {x k } k∈N is a sequence generated by the SKM method and Assumptions 1 and 2 hold and {µ k } k∈N is generated using a scheme outlined in Table 2 , then Assumption 3 holds and, by Theorem 3.1, the limit (33) holds.

4. A NEURAL NETWORK VIEW:

The SKM method may be truncated and executed inferences with a neural network. The input into the network is the data d, often in vector form. Each layer of the network consists of an iterate x k , with the output from the final layer K providing the inference estimate x K that approximates a solution to (1). The feedforward operations between layers are formed by the operator T L2O (•, ζ k ) if the descent condition in Line 6 of the SKM method holds and T otherwise. We encode all the network parameters with Θ := (ζ k ) K k=1 . The set C, over which Θ is minimized, may be chosen with great flexibility, the only restriction being that C is a subset of the K-tuple set where each ζ k is such that T L2O (•; ζ k ) is well-defined. Since C determines the form of possible Θ, we refer to it in Algorithm 2 as the network structure. For each application of the algorithm, the operators T and S change, depending upon the particular data d. For example, the data d could correspond to the measurement vector d = Ax when attempting to recover x. Thus, to make explicit this dependence of each operator on the data d, we henceforth incorporate a subscript to write T d and S d .

The "optimal" choice of parameters Θ depends upon the particular application. Suppose each d is drawn from a common distribution D. Then a choice of "optimal" parameters Θ may be identified as those for which the expected value of φ d (x K ) is minimized, where φ d : H → R is an appropriate cost function. In mathematical terms, this is expressed by stating Θ forms a solution to the problem

EQUATION

where the expectation E d∼D is taken over all d ∈ D and we emphasize the dependence of x K on Θ and d by writing

x K = x K (Θ, d).

In practice, the problem ( 13) is approximately solved by using a sample set of data {d n } N n=1 taken from the distribution D and minimizing an empirical loss function. We summarize the procedure for creating, training, and performing inferences with such L2O networks in Algorithm 2. There M denotes the neural network, which is dependent upon Θ, δ, and the choice of safeguarding used to construct each µ k (see Line 6). Remark 4.1. Different learning problems than (13) may be used (e.g., the min-max problem used by adversarial networks (Goodfellow et al., 2014) ). Remark 4.2. Two particular choices of the set C are of interest. The first case is when each parameter ζ k may be chosen independently, i.e., the network weights vary by layer. This is used in our numerical examples below. The second case is when C is restricted such that each of its K entries are identical, i.e., the parameters across all layers are fixed so that

ζ 1 = ζ 1 = • • • = ζ K .

This latter case corresponds to the structure of recurrent neural networks (RNNs).

5. NUMERICAL EXAMPLES:

This section presents example L2O schemes applied within the LSKM framework. 1 We numerically investigate (i) the convergence rate of LSKM relative to corresponding KM iterations, (ii) the efficacy of safeguarding procedures when inferences are performed on data for which L2O fails intermittently, and (iii) the convergence of LSKM schemes even when the application of T L2O is not justified theoretically. We first use T L2O from ALISTA (Liu et al., 2019a) on the LASSO Problem. Then we apply the differentiable linearized ADMM of Xie et al. (2019) to a similar problem. We also implement a new L2O scheme for solving the nonnegative least squares (NNLS) problem.

In all experiments, we take µ 1 = S d (x 1 ) 2 . Implementations of the LSKM Algorithm are abbreviated by the scheme in Table 2 , e.g., we write GS(α) to mean the Geometric Sequence method with parameter α is used to construct {µ k } k∈N . As in Algorithm 2, training is completed without safeguarding. For notational compactness, in each example we let f d the optimal value of f d (x) among all possible x. Note each objective function is dependent upon the data d, which implies it it is only meaningful to compare curves below that are within the same plot (and not appropriate to compare curves in separate plots). We illustrate the performance in plots using an approximation of the expected objective error

EQUATION

using 1,000 test samples. The optimal value f d for each sample d is estimated by running the KM method for 15,000 iterations.

5.1. ALISTA FOR LASSO:

Here we consider the popular LASSO problem for sparse coding. Consider an unknown sparse vector x ∈ R n and a matrix A ∈ R m×n , which is called the dictionary. Then assume we have access to noisy linear measurements d ∈ R m , where ε ∈ R m is additive Gaussian white noise and

EQUATION

Even for underdetermined systems, when x is sufficiently sparse and τ ∈ (0, ∞) is an appropriately chosen regularization parameter, x can often be recovered faithfully by solving the LASSO problem

EQUATION

where • 2 and • 1 are the 2 and 1 norms, respectively. A classic method for solving (16) is the iterative shrinkage thresholding algorithm (ISTA) (e.g., see (Daubechies et al., 2004 )), a special case of the proximal-gradient method in Table 1 . Given x 1 ∈ R n , this method iteratively computes

EQUATION

where L = A t A 2 and η θ is the soft-thresholding function defined by component-wise operations:

EQUATION

)

Implementation details for the LASSO problem may be found in the Appendix.

We applied the LSKM Algorithm to the LASSO problem above by using T defined in ( 17) and the update operation from ALISTA (Liu et al., 2019a) . The tunable operator T L2O is parameterized by ζ = (θ, γ) for positive scalars θ and γ and defined by

EQUATION

with W defined in the Appendix. The parameter Θ used in Algorithm 2 is Θ = (θ k , γ k ) K k=1 , which consists of 2K scalars. Note T L2O may fail to be averaged, depending on the choice of ζ = (θ, γ).

The primary illustration of the rapid convergence using LSKM-ALISTA relative to the KM counterpart ISTA is in Figure 1a . There each x k estimating a solution to ( 16) is computed for data d drawn from the same distribution D s that was used to train the LSKM network. Figure 1b shows a plot comparing performance of LSKM and KM with each d there drawn from a distribution D u that is different than D s . For this reason, we refer to D u as the unseen distribution. For d ∈ D u , the L2O scheme ALISTA usually fails to converge without safeguarding whereas the safeguarded method LSKM-ALISTA still converges and much faster than ISTA. The dotted plot with square markers shows the percentage of safeguard activations that occurred in texting. The only two layers in Figure 1b where safeguarding kicked in were layers 2, 7 and 13, with respective percentages 36.7%, 99.6% and 99.2% among 1000 testing samples. Additionally summary results are provided in Table 3 . These show the relative usefulness of difference choices of training loss function φ d , safeguarding method, and performance measured by different test loss functions. In particular, using φ d = f d and EMA(0.25) safeguarding yielded the lowest function value f d ; however, instead using GS(0.1) safeguarding resulted in estimates closer to the underlying sparse signals x . In both the seen and unseen cases, 20 iterations of LSKM-ALISTA yields better function values f d than thousands of iterations of ISTA, which reveals orders of magnitude speedup by the safeguarded L2O schemes.

5.2. LINEARIZED ADMM:

Let A ∈ R m×n , x ∈ R n , and d ∈ R m be as in Subsection 5.1. Here we apply the L2O scheme differentiatable linearized ADMM (D-LADMM) of Xie et al. (2019) to the closely related sparse coding problem

EQUATION

)

The operators defining T and T L2O for LSKM-D-LADMM are provided in the Appendix along with further implementation details. Comparison plots are provided in Figure 2 and Table 4 summarizes further results.

5.3. PROJECTED GRADIENT FOR NONNEGATIVE LEAST SQUARES:

Let A ∈ R m×n and d ∈ R m . Here we consider the nonnegative least squares (NNLS) problem

EQUATION

)

We proceed by using the projected gradient method in Table 1 , where

C := {x ∈ R n : x ≥ 0}, ∇f d (x) = A T (Ax -d)

, and proj C (x) = max(x, 0), with the max applied component-wise. Then, for α = 1/ A T A 2 , we take

EQUATION

)

10 0 10 1 10 2 10 3 10 4 10 -5 We then relax T to obtain, for ζ = (α, β, W ) with α ∈ R, W ∈ R n×n , and β ∈ R n , to obtain

10 -2 10 1 Iteration/Layer k E d∼Ds [f d (x k ) -f d ] ISTA LSKM-ALISTA (a)

EQUATION

In the special case that W = A T , α = 1/ A T A 2 , and β = 0, we recover ( 22). When safeguarding does not kick in, ( 23) yields a familiar network structure with the max as the activation function.

Here

Θ = (ζ k ) K k=1 = (W k , D k , β k ) K k=1 consists of (n 2 + 2n)K trainable parameters.

During training we force β K ≥ 0 to ensure x K ∈ C. However, note x k might not be in C for k < K. In this example, each learned W k is likely an approximation of the pseudo-inverse of A. Summary plots are given in Figure 3 and more results in Table 5 of the Appendix.

6. DISCUSSION:

Our numerical examples provide several insights. Tables 3 to 5 in the Appendix reveal the choice of training loss function φ d results in different models and performance, even though the minimizers of each φ d are identical. This is due to the difference in gradient directions of each φ d and reveals the 'best' choice of φ d may depend on the application. Figure 1b shows the necessity of safeguarding to ensure convergence and that inferences on unseen data can still be much more rapid than traditional KM methods. For the space limit, we do not focus on timed results. We simply note that the relative cost per iteration of the SKM method versus KM is roughly double in the worst case (since two tentative update operations are needed). However, in cases like our ALISTA example, this relative cost is greatly reduced because some computations (e.g., Ax k -d) are common to both T L2O and T . In each example, L2O schemes obtain an order of magnitude reduction in computational cost.

The results of this work reveal a practical procedure for creating and implementing L2O schemes on new problems. Suppose one is presented a convex optimization problem and a standard general purpose method for solving it. The method determines an averaged operator T (e.g., as in Table 1). To construct T L2O , one may let each scalar, vector, and matrix that does not represent the data d be parameterized (i.e., all of its entries learnable). We may further the operations by replacing scalar parameters with vectors and using element-wise products as in D-LADMM (see (31)). Having T and T L2O , a network can be constructed via the LSKM algorithm. For training the network, the learnable parameters may be initialized to the quantities that reduce T L2O to the original operator T . Additional structures of parameters may be invoked (e.g., as done by ALISTA with W in ( 26)). The number of terms to parameterize and how to structure parameterizations are matters subject to the a practitioner's needs/priorities for a given application.

10 0 10 1 10 2 10 3 10 -1 10 0 10 1 10 2 Iteration/Layer k E d∼Ds [f d (x k ) -f d ]

7. CONCLUSION:

This work establishes a framework for extending the L2O methodology to a wide class of iterative optimization procedures (i.e., KM methods). We provide theoretical results that demonstrate sequences generated by our SKM method possess the property that the distance between each iterate and fixed point set converges to zero. When there is a unique cluster point, this is equivalent to convergence of the method to a solution. Furthermore, our LSKM algorithm provides a straightforward neural network design that inherits the theoretical properties from the SKM method. Practical guidance is also given for constructing learnable operations for new problems. Our numerical experiments demonstrate order(s) of magnitude faster convergence by LSKM implementations than general-purpose counterparts and the efficacy of safeguarding in cases where L2O schemes would otherwise fail to converge. Future work will provide a more efficacious fall-back method for using data-driven updates than classic KM updates and will investigate stochastic extensions with quasi-Féjer monotone operators. A.1 NUMERICAL EXAMPLE SUPPLEMENT MATERIALS Supplement to Subsection 5.1. In similar manner to (Chen et al., 2018) and (Liu et al., 2019a) , we use the following set up. We take m = 250, n = 500, and τ = 0.001. Each entry of the dictionary A is sampled i.i.d from the standard Gaussian distribution, i.e., a ij ∼ N (0, 1/m). Having these entries, we then normalize each column of A, with respect to the Euclidean norm. Each d in the distribution D s of data used to train the neural network is constructed using (15) with noise ε ∼ 0.1•N (0, 1/m) and each entry of x as the composition of Bernoulli and Gaussian distributions, i.e., x j ∼ Ber(0.1) • N (0, 1) for all j ∈ [n]. Each d in the unseen distribution D u is computed using the same distribution of noise ε as before and using x j ∼ Ber(0.2) • N (1, 2). Our data set consists of 10,000 training samples and 1,000 test samples.

If τ is chosen too large, then the solution to the problem ( 16) is simply the zero vector, which is undesirable. If τ is chosen too small, then the solution will not be sparse. Thus, appropriate choice of τ is crucial. Since the goal of the LASSO problem is to faithfully recover the underlying sparse vector x used to create the measurement d, we find it reasonable to choose τ as a solution estimate of the problem min τ ∈(0,∞)

EQUATION

where x K d is the output of ISTA after K = 15, 000 iterations for the problem ( 16) and x d is the sparse vector from which d was created. To simplify (24), we assume τ was of the form t = 10 α with integer α ∈ Z to obtain the problem

EQUATION

Approximately solving this discrete problem for our data set revealed the optimal choice yielded α = -3 and τ = 0.001.

In practice, such a choice of τ may be difficult to ascertain. However, if the practitioner has access to each underlying x d for their training data d, then one can circumvent the problem of choosing τ by simply using the training loss function φ

EQUATION

2 . This was precisely the approach taken in Liu et al. (2019a) .

To define the learned operations in T L2O , we let

EQUATION

where • F is the Frobenius norm and the Matlab notation M :, is used to denote the th column of the matrix M .

Supplement to Subsection 5.2. LADMM is used to solve problems of the form

EQUATION

for which LADMM generates sequences {x k } k∈N , {z k } k∈N and {ν k } k∈N defined by the updates 

EQUATION

NA f d f d f d f d R f,Ds (x K ) 3.48e+00 1.78e+00 1.52e+01 1.78e+00 1.28e+01

with given scalars α, β, γ ∈ (0, ∞). The problem (20) may be written in the form of ( 27) by taking 

f = τ • 1 , g = • 1 ,

NA f d f d f d f d R f,Ds (x K ) 3.80e+06 1.86e+06 2.37e+07 1.86e+06 1.65e+08

And, if {x k } k∈N contains a single cluster point, then {x k } k∈N converges to a point x ∈ Fix(T ).

Proof. If the inequality in Line 6 holds finitely many times, then there exists an index beyond which Line 9 is always used to update x k . In this case, for large k the SKM Method takes the form of the classic KM Method, which is known to converge (e.g., see (Cegielski, 2012, Theorem 3.7 .1), (Groetsch, 1972, Corollary 3) , and (Bauschke & Combettes, 2017, Theorem 5.15) ). Thus, it is sufficient to only consider the case where the inequality in Line 6 holds infinitely many times. We proceed by first showing the sequence {x k } k∈N is bounded (Step 1). This is used to prove {x k } k∈N has a cluster point in Fix(T ) (Step 2). Results from these steps are then applied to obtain the desired limit (33) (Step 3).

Step 1: By Assumption 2, there exists R ∈ (0, ∞) sufficiently large to ensure

x > R =⇒ S(x) > 1, for all x ∈ H.

Equivalently, we may write S(x) ≤ 1 =⇒ x ≤ R, for all x ∈ H.

By Assumption 3, there also exists N 1 ∈ N such that µ k ≤ 1, for all k ≥ N 1 .

(36)

Fix any z ∈ Fix(T ). We claim

x k -z ≤ max ∈[N1]

{2R, x -z }, for all k ∈ N.

(37)

The result (37) holds trivially for all k ∈ [N 1 ]. Proceeding by induction, suppose (37) holds for some k > N 1 . If the inequality in Line 6 holds, then (35), ( 36) and the update formula in Line 7 together imply x k+1 ≤ R. Since S(z) = 0, (35) also implies z ≤ R. Thus,

EQUATION

)

If instead the update in Line 9 is applied, the averagedness of T implies there is α ∈ (0, 1) such that

EQUATION

(e.g., see Prop. 4.35 in (Bauschke & Combettes, 2017) or Cor. 2.2.15 and Cor. 2.2.17 in (Cegielski, 2012) ), and so

x k+1 -z ≤ x k -z ≤ max ∈[N1]

{2R, x -z }.

(40)

Equations ( 38) and ( 40) together close the induction, from which (37) follows. Whence

x k ≤ x k -z + z ≤ max ∈[N1]

{2R, x -z } + R, for all k ∈ N, (41) which verifies the sequence {x k } k∈N is bounded.

Step 2: Because the inequality in Line 6 holds infinitely many times, there exists a subsequence {x q k } k∈N ⊆ {x k } k∈N satisfying

EQUATION

from which the squeeze theorem asserts S(x q k ) → 0. Since {x k } k∈N is bounded, so also is {x q k } k∈N . Thus, there exists a subsequence {x k } k∈N ⊆ {x q k } k∈N converging to a limit p ∈ H. Then applying the fact S is 2-Lipschitz and • is continuous yields

EQUATION

That is, {x k } k∈N contains a cluster point p ∈ Fix(T ).

Step 3: Let ε > 0 be given. Following (Bauschke et al., 2014, Def. 2), define the ε-enlargement 

where the first implication holds because x ∈ B(0, R) by ( 35) and x / ∈ C by (46). By Assumption 3, there exists N 2 ∈ N such that

EQUATION

)

By the result of Step 2, there exists N 3 ≥ N 2 such that

EQUATION

)

We claim d Fix(T ) (x k ) ≤ ε, for all k ≥ N 3 , (50) which, by the arbitrariness of ε, implies (33) holds. Indeed, inductively suppose (50) holds for some k ≥ N 3 . If the inequality in Line 6 holds, then (47) and (48) together imply

EQUATION

Otherwise, letting P : H → H be the projection operator onto Fix(T ), we deduce

EQUATION

:

All of the codes for this work can be found on Github here: (link will be added after review process).

Figures and Tables:

Learned SKM (LSKM) 1: Stage 1: Initialization/Training. 2: Choose envelope T L2O (•; •) and network structure C, parameterized by Θ = (ζ k ) K k=1 3: Choose training loss function φ d 4: Choose 'optimal' parameter Θ ∈ arg min Θ∈C E d∼D φ d (x K ) , assuming µ k = ∞ at each layer k 5: Choose δ and safeguarding scheme for {µ k } K k=1 6: Define the neural network M = M Θ ,δ,µ k . 7: Stage 2: Inference. 8: For input d return x = M(d)

Performance on seen distribution, i.e., d ∼ Ds Performance on unseen distribution, i.e., d ∼ Du

Figure 1: Expected function value error versus iteration when applied to two different data distributions. Training used φ d = f d . Inferences used δ = 0.01 and EMA(0.1).

Performance on unseen distribution, i.e., d ∼ Du

Figure 2: Expected relative function value error versus iteration when applied to two different distributions. Training used φ d = f d . Inferences used δ = 0.01, and GS(0.1) in (a) and EMA(0.75) in (b).

Figure 3: Expected function value error versus iteration when applied to two different data distributions. Training used φ d = f d . Inferences used δ = 0.01 and GS(0.1).

(T ) [ε] := {x ∈ H : d Fix(T ) (x) ≤ ε}. (44) Note Fix(T ) [ε] is a nonempty closed and bounded subset of H. SetC := B(0, R) -Fix(T ) [ε/2] ∪ ∂Fix(T ) [ε/2] ,(45)where B(0, R) is the closed ball of radius R centered at the origin. Because H is finite dimensional and C closed and bounded, C is compact. Thus, every continuous function obtains its infimum over C. In particular, we may setζ = min x∈C S(x) . (46) Note ζ > 0 since C ∩ Fix(T ) = ∅. Consequently, letting ζ := min{1, ζ/2} yields S(x) ≤ ζ =⇒ x ∈ Fix(T ) [ε] =⇒ d Fix(T ) (x) ≤ ε, for all x ∈ H,

Summary of LASSO Problem Results on data d drawn from the seen distribution D s . K = 10, 000 iterations of ISTA are used by the KM reference method 'ISTA 10K' while K = 20 and δ = 0.01 for the LSKM schemes. 'No guard' refers to LSKM without safeguarding.

Table content: <table><tr><td>Method</td><td colspan="7">ISTA 10K No guard EMA(0.1) EMA(0.1) EMA(0.1) GS(0.1)</td><td>RM(3)</td></tr><tr><td>Loss φ d</td><td>NA</td><td>f d</td><td>1 2</td><td>S d (•) 2 2</td><td>S1(•) 1</td><td>f d</td><td>f d</td><td>f d</td></tr><tr><td>R f,Ds (x K )</td><td>7.72e-04</td><td>3.36e-04</td><td colspan="2">4.65e-04</td><td>4.58e+00</td><td>3.33e-04</td><td colspan="2">3.32e-04 6.68e-04</td></tr><tr><td>E[ 1 2 S(x K ) 2 2 ]</td><td>2.11e-09</td><td>4.68e-07</td><td colspan="2">3.75e-07</td><td>1.58e-05</td><td>4.68e-07</td><td colspan="2">4.68e-07 4.95e-07</td></tr></table>

Summary of D-LADMM Problem Results on data d drawn from the seen distribution D s . K = 1, 000 iterations of LADMM are used by the KM reference method 'LADMM 1K' while K = 20 and δ = 0.01 for the LSKM schemes. 'No guard' refers to LSKM without safeguarding.

Table content: <table><tr><td>Method</td><td>LADMM 1K No guard EMA(0.75)</td><td>GS(0.1)</td><td>RT</td></tr><tr><td>Loss φ d</td><td/><td/><td/></tr></table>

and B = -Id. In this case, the proximal operators in (28) reduce to softthresholding operators. Although not given in Table1, the update ν k+1 is generated by applying an averaged operator T to ν k . (This follows since LADMM may expressed as a special case of proximal ADMM, which itself is a special case of the ADMM method. The details of this derivation are outside the scope of this paper.) This implies 10 0 10 1 10 2 10 3 b) Performance on unseen distribution, i.e., d ∼ Du

Table content: <table><tr><td>E d∼Ds [f d (x k ) -f d ]</td><td>10 0 10 2 10 4</td><td>PG LSKM-PG</td><td>E d∼Du [f d (x k ) -f d ]</td><td>10 0 10 2 10 4</td><td/><td/><td>PG LSKM-PG L2O-PG</td><td>20% 40% 60% 80%</td><td>Safeguard Activation Freq.</td></tr><tr><td/><td>10 -2</td><td/><td/><td>10 0 10 -2</td><td>10 1</td><td>10 2</td><td colspan="2">10 3 0%</td></tr><tr><td/><td/><td>Iteration/Layer k</td><td/><td/><td/><td>k</td><td/></tr><tr><td colspan="3">(a) Performance on seen distribution Ds</td><td/><td/><td/><td/><td/></tr></table>

Summary of NNLS Problem Results on data d drawn from the seen distribution D s . K = 1, 500 iterations of PG are used by the KM reference method 'PG 1.5K' while K = 20 and δ = 0.01 for the LSKM schemes. 'No guard' refers to LSKM without safeguarding.

Table content: <table><tr><td>Method</td><td>PG 1.5K No guard EMA(0.75)</td><td>GS(0.1)</td><td>RT</td></tr><tr><td>Loss φ d</td><td/><td/><td/></tr></table>

Peer review (iclr2020):
This paper presents a unified framework for parametrizing provably convergent algorithms and learning the parameters for a training dataset of problem instances of interest. The learned algorithm can then be used on unseen problems. One key idea to this algorithm is that it is safeguarded, meaning it will perform some standard, non-learned iterations, if the predicted iterate is not good enough under some condition. There are three main features of the proposed approach: 1- It unifies various previous approaches such as LISTA, ADMM, non-negative least squares, etc. By defining some operators and safeguarding rules, the same learning approach can be leveraged for these different optimization problems. 2- It is shown that the learned algorithms are provably convergent under some mild assumptions. 3- Empirically, it is shown that the learned algorithms converge faster than the non-learned counterpart on sparse coding, ADMM and non-negative least squares; they use safeguarding sparingly, particularly when used to solve test instances from the same distribution as the training instances. Additionally, the paper is very well-written. I did not verify the proofs in detail but they seemed OK at a high-level; however, I am not an expert in convex optimization so I hope other reviewers will be able to comment on this aspect. I do have some deep concerns about the evaluation metrics used to report the results that I will discuss next; these are the main reason for my current score, but I am willing to adjust it if the authors address them convincingly. I also have some comments about related work. Experimental evaluation: - The error metric (15) is not suitable for evaluating the performance of an optimization algorithm. You should compute the expectation of the relative error, i.e. E_{d~D} [(f_d(x)-f_d^*) / f_d^*]. This is similar to the average approximation ratio used in the learning to optimize papers for discrete problems (see refs. below). (15) is just the ratio of the expected absolute error to the expected optimal value; I don*t think that is equivalent to what I suggested. - The relative error values are massive in some cases, e.g. Fig. 3. What*s going on there? Are all methods performing that horribly? Am I misinterpreting the metric? - Why do the plots for the seen distribution extend over thousands of iterations but only for tens of iterations for the unseen distribution? - Please use the same scale for the y-axes in Figs. 1-3. Methodology: - Your method requires learning per-iteration parameters. The other L2O methods for gradient descent (see refs. below) use shared parameters instead. This allows them to run for many iterations, possibly beyond what they were trained for. Your method does not allow for that. On the other hand, such models are recurrent and thus possibly more difficult to train than your unrolled feedforward model. Is the fixed number of iterations a limitation of your method? Please discuss this. Related work: - Learning for gradient descent: I am surprised these papers are not mentioned although they are quite relevant. They are rather recurrent networks with shared parameters across iterations, but you should also compare against them both conceptually and experimentally: *Learning to optimize.* arXiv preprint arXiv:1606.01885 (2016). *Learning to learn by gradient descent by gradient descent.* Advances in neural information processing systems. 2016. - Learning to optimize in the discrete setting: there is lots of recent work on this that you should at least point to in passing, e.g.: *Learning combinatorial optimization algorithms over graphs.* Advances in Neural Information Processing Systems. 2017. *Combinatorial optimization with graph convolutional networks and guided tree search.* Advances in Neural Information Processing Systems. 2018. (Survey) *Machine Learning for Combinatorial Optimization: a Methodological Tour d*Horizon.* arXiv preprint arXiv:1811.06128 (2018). - Theory for learning to optimize: Since you have a theoretical basis for your framework, you should discuss connections to other recent frameworks such as the one below by Balcan et al. It is geared towards the discrete setting and sample complexity rather than convergence, but you should nevertheless discuss it. Balcan, Maria-Florina, et al. *How much data is sufficient to learn high-performing algorithms?.* arXiv preprint arXiv:1908.02894 (2019). Clarification questions: - *The choice of parameter ? k in Line 3 may be any value that results in a well-de?ned operator T L2O*: what is *well-defined* here? that T_{L20} is averaged? Minor: - Page 3: *A classic theorem states sequences* -> *A classic theorem states that sequences* - Appendix proofs: please organize into sections and restate the statements before the proofs.