SAFEGUARDED LEARNED CONVEX OPTIMIZATION
Year: 2019

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 DLADMM) 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.

heading: 1 INTRODUCTION
text: 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.

heading: 2 FIXED POINT METHODS
text: Let H be a finite dimensional Hilbert space (e.g., the Euclidean space Rn) with inner product 〈·, ·〉 and norm ‖ · ‖. Denote the set of fixed points of each operator T : H → H by Fix(T ) := {x ∈ H : Tx = x}. In this work, for an operator T with a nonempty fixed point set (i.e., Fix(T ) 6= ∅), 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 : H → H is nonexpansive if ‖Tx− Ty‖ ≤ ‖x− y‖, for all x, y ∈ H. (2)
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
dC(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
Rα∂f := 2Jα∂f − Id. (5) If f is closed, convex, and proper, then the resolvent is precisely the proximal operator, i.e.,
Jα∂f = proxαf (x) := arg min z∈H
αf(z) + 1
2 ‖z − x‖2. (6)
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 {xk}k∈N with arbitrary x1 ∈ H satisfies the update relation
xk+1 = T (xk), for all k ∈ N, (7) then there is x? ∈ Fix(T ) such that {xk}k∈N converges to x?, i.e., xk → x?.
There are pathological cases where the result fails for operators that are only nonexpansive (e.g., when x1 6= 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.

heading: 3 SAFEGUARDED KM METHOD
text: This section generalizes the classic KM iteration in (7). We accomplish this by defining an envelope of operators TL2O(·; ·). For a parameter ζ chosen from an appropriate set, we let TL2O(·; ζ) define an operator on H. Changing ζ may define a new operator with different properties. We do not impose restrictions on TL2O(·; ζ) other than it be well-defined, meaning TL2O(·; ζ) may fail to be averaged and/or fail to have a fixed point. This is illustrated by the following two examples.
Example 3.1. Let Q : H → H be nonexpansive. Then define TL2O : H× R→ R by TL2O(x; ζ) := (1− ζ)x+ ζQ(x). (8)
For ζ ∈ (0, 1), the operator TL2O(·; ζ) 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)). 4
Example 3.2. Let f : H → R be closed, convex and proper, and define TL2O : H× (0,∞)→ R by
TL2O(x; ζ) := proxζf (x). (9)
For fixed ζ ∈ (0,∞), the operator TL2O(·; ζ) is averaged. 4
The practicality of TL2O 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 x1 = 0. The for loop on Lines 2 to 11
generates each update xk+1 in the sequence {xk}k∈N. The choice of parameter ζk in Line 3 may be any value that results in a well-defined operator TL2O(·; ζ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, TL2O(xk; ζk) is used to update xk via Line 7. Otherwise, a KM update is used to update xk 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.,
lim ‖x‖→∞
‖S(x)‖ =∞. (10)
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.
Algorithm 2 Learned SKM (LSKM) 1: Stage 1: Initialization/Training. 2: Choose envelope TL2O(·; ·) and network structure C, parameterized by Θ = (ζk)Kk=1 3: Choose training loss function φd 4: Choose ‘optimal’ parameter
Θ? ∈ arg min Θ∈C
Ed∼D [ φd(x K) ] ,
assuming µk =∞ at each layer k 5: Choose δ and safeguarding scheme for {µk}Kk=1 6: Define the neural networkM =MΘ?,δ,µk . 7: Stage 2: Inference. 8: For input d return x =M(d)
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 {xk}k∈N is a sequence generated by the SKM method and Assumptions 1 to 3 hold, then
lim k→∞
dFix(T ) ( xk ) = 0. (11)
And, if {xk}k∈N contains a single cluster point, then {xk}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 xk and (possibly) previous iterates. These update schemes enable each µk to trail the value of the residual norm ‖S(xk)‖. This implies there may exist j ∈ N for which
‖S(xj)‖ < ‖S(xj+1)‖ = ‖S ( TL2O(x j ; ζj) ) ‖ ≤ (1− δ)µj . (12)
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 {xk}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.

heading: 4 A NEURAL NETWORK VIEW
text: 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 xk, with the output from the final layer K providing the inference estimate xK that approximates a solution to (1). The feedforward operations between layers are formed by the operator TL2O(·, ζk) if the descent condition in Line 6 of the SKM method holds and T otherwise. We encode all the network parameters with Θ := (ζk)Kk=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 TL2O(·; ζ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 Td and Sd.
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(xK) 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
min Θ
Ed∼D [ φd(x K(Θ, d)) ] , (13)
where the expectation Ed∼D is taken over all d ∈ D and we emphasize the dependence of xK on Θ and d by writing xK = xK(Θ, d). In practice, the problem (13) is approximately solved by using a sample set of data {dn}Nn=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. ThereM 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).

heading: 5 NUMERICAL EXAMPLES
text: 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 TL2O is not justified theoretically. We first use TL2O 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 = ‖Sd(x1)‖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 fd(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
Ed∼D [ fd(x k)− f?d ] , (14)
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.

heading: 5.1 ALISTA FOR LASSO
text: Here we consider the popular LASSO problem for sparse coding. Consider an unknown sparse vector x? ∈ Rn and a matrix A ∈ Rm×n, which is called the dictionary. Then assume we have access to noisy linear measurements d ∈ Rm, where ε ∈ Rm is additive Gaussian white noise and
d = Ax? + ε. (15)
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
min x∈Rn
fd(x) := 1
2 ‖Ax− d‖22 + τ‖x‖1, (16)
1All of the codes for this work can be found on Github here: (link will be added after review process).
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 x1 ∈ Rn, this method iteratively computes
xk+1 := T (xk) := ητ/L ( xk − 1
L AT (Axk − d)
) , for all k ∈ N, (17)
where L = ‖AtA‖2 and ηθ is the soft-thresholding function defined by component-wise operations:
ηθ(x) := sgn(x) ·max{0, |x| − θ}. (18)
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 TL2O is parameterized by ζ = (θ, γ) for positive scalars θ and γ and defined by
TL2O(x; ζ) := ηθ ( x− γWT (Ax− d) ) , (19)
with W defined in the Appendix. The parameter Θ used in Algorithm 2 is Θ = (θk, γk)Kk=1, which consists of 2K scalars. Note TL2O 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 xk estimating a solution to (16) is computed for data d drawn from the same distribution Ds 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 Du that is different than Ds. For this reason, we refer to Du as the unseen distribution. For d ∈ Du, 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 = fd and EMA(0.25) safeguarding yielded the lowest function value fd; 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 fd than thousands of iterations of ISTA, which reveals orders of magnitude speedup by the safeguarded L2O schemes.

heading: 5.2 LINEARIZED ADMM
text: Let A ∈ Rm×n, x? ∈ Rn, and d ∈ Rm 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
min x∈Rn
‖Ax− d‖1 + τ‖x‖1. (20)
The operators defining T and TL2O 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.

heading: 5.3 PROJECTED GRADIENT FOR NONNEGATIVE LEAST SQUARES
text: Let A ∈ Rm×n and d ∈ Rm. Here we consider the nonnegative least squares (NNLS) problem
min x∈Rn
fd(x) := 1
2 ‖Ax− d‖22 s.t. x ≥ 0. (21)
We proceed by using the projected gradient method in Table 1, where C := {x ∈ Rn : x ≥ 0}, ∇fd(x) = AT (Ax− d), and projC(x) = max(x, 0), with the max applied component-wise. Then, for α = 1/‖ATA‖2, we take
T (x) := max ( x− αAT (Ax− d), 0 ) . (22)
We then relax T to obtain, for ζ = (α, β,W ) with α ∈ R, W ∈ Rn×n, and β ∈ Rn, to obtain TL2O(x; ζ) := max ( x− αWT (Ax− d), β ) . (23)
In the special case that W = AT , α = 1/‖ATA‖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)Kk=1 = (W k, Dk, βk)Kk=1 consists of (n 2 + 2n)K trainable parameters. During training we force βK ≥ 0 to ensure xK ∈ C. However, note xk 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.

heading: 6 DISCUSSION
text: 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., Axk − d) are common to both TL2O 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 TL2O, 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 generalize the operations by replacing scalar parameters with vectors and using element-wise products as in D-LADMM (see (31)). Having T and TL2O, a network can be constructed via the LSKM algorithm. For training the network, the learnable parameters may be initialized to the quantities that reduce TL2O 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.

heading: 7 CONCLUSION
text: 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 quasiFéjer monotone operators.

heading: A APPENDIX
text: This Appendix contains a subsection with supplemental materials and a subsection with proofs of the theoretical results.
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 dictionaryA is sampled i.i.d from the standard Gaussian distribution, i.e., aij ∼ 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 Ds 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 distributionDu 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,∞)
Ed∼D [ ‖xKd − x?d‖22 ] , (24)
where xKd 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
min α∈Z
Ed∼D [ ‖xKd − x?d‖22 ] , (25)
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 φd(x) := ‖x − x?d‖22. This was precisely the approach taken in Liu et al. (2019a).
To define the learned operations in TL2O, we let
W ∈ arg min M∈Rm×n ‖MTA‖F , s.t.(M:,`)TA:,` = 1, for all ` ∈ [n], (26)
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
min x∈Rn,z∈Rm f(x) + g(z) s.t. Ax+Bz = d, (27)
for which LADMM generates sequences {xk}k∈N, {zk}k∈N and {νk}k∈N defined by the updates xk+1 := proxβf ( xk − βAT [ νk + α ( Axk +Bzk − d )]) ,
zk+1 := proxγg ( zk − γBT [ νk + α ( Axk+1 +Bzk − d )]) ,
νk+1 := νk + α ( Axk+1 +Bzk+1 − d ) ,
(28)
with given scalars α, β, γ ∈ (0,∞). The problem (20) may be written in the form of (27) by taking f = τ‖ · ‖1, g = ‖ · ‖1, and B = −Id. In this case, the proximal operators in (28) reduce to softthresholding operators. Although not given in Table 1, 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
‖S(νk)‖ = ‖νk − T (νk)‖ = ‖νk − νk+1‖ = α‖Axk+1 − zk+1 − d‖. (29) Because we compare methods that use differing choices of α, the residual comparison illustrations presented below will, for consistency of interpretation, assume α = 1 in (29). And, although xk+1 and zk+1 form intermediate computations, for notational clarity, the term xk in the SKM and LSKM schemes is replaced in this subsection by the tuple (xk, zk, νk). This is of practical importance too since it is the sequence {xk}k∈N that converges to a solution of (20). We now modify the iteration (28) for the problem (20) to create the D-LADMM L2O scheme. We generalize soft-thresholding to vectorized soft-thresholding for β ∈ Rn by
ηβ(x) = (ηβ1(x1), ηβ2(x2), . . . , ηβn(xn)). (30)
We assume ηβ represents the scalar soft-thresholding in (18) when β ∈ R and the vector generalization (30) when β ∈ Rn. Combining ideas from ALISTA (Liu et al., 2019a) and D-LADMM (Liu et al., 2019b), given (xk, zk, νk) ∈ Rn × Rm × Rm, αk, γk, ξk ∈ Rm, βk, σk ∈ Rn, W1 ∈ Rn×m, and W2 ∈ Rm×m, set
x̃k+1 := ηβk ( xk − σk ◦ (W k1 )T [ νk + αk ◦ ( Axk − zk − d )]) ,
z̃k+1 := ηγk ( zk − ξk ◦ (W k2 )T [ νk + αk ◦ ( Ax̃k+1 − zk − d )]) ,
ν̃k+1 := νk + αk ◦ ( Ax̃k+1 − z̃k+1 − d ) ,
(31)
with element-wise products denoted by ◦. For the parameter ζk := (αk, βk, γk, σk, ξk,W k1 ,W k2 ), then define
TL2O(x k, zk, νk; ζk) := (x̃k+1, z̃k+1, ν̃k+1). (32)
Fixing the number of iterations K, the learnable parameters from (31) used in the LSKM Algorithm may be encoded by Θ = (ζk)Kk=1 = ( αk, βk, γk, σk, ξk,W k1 ,W k 2 )K k=1
, consisting of (2n+ 3m+ mn+m2)K scalars.
For data generation, we follow the same settings as in the experiments of Subsection 5.1. Although the value of τ is not optimally chosen, this was used as it worked well for D-LADMM.
Supplement to Subsection 5.3. We take m = 500, n = 250, aij ∼ Ber(0.1) · rand[0, 1]. We use noise ε ∼ 0.1 · N (0, 1/m). Each d ∼ Ds used for training the neural network is sampled using d = Ax? + ε with x ∼ max(N (0, 1), 0). For unseen data d ∼ Du we sample x? ∼ max(N (5, 102), 0). We sample 10,000 training samples from Ds to train the neural network and 1,000 samples from Ds and Du, respectively, for testing.
A.2 PROOFS
For ease of reference, we restate Theorem 3.1 below. Then provide a proof is provided for it.
Theorem 3.1. If {xk}k∈N is a sequence generated by the SKM method and Assumptions 1 to 3 hold, then
lim k→∞
dFix(T ) ( xk ) = 0. (33)
And, if {xk}k∈N contains a single cluster point, then {xk}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 xk. 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 {xk}k∈N is bounded (Step 1). This is used to prove {xk}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. (34)
Equivalently, we may write
‖S(x)‖ ≤ 1 =⇒ ‖x‖ ≤ R, for all x ∈ H. (35) By Assumption 3, there also exists N1 ∈ N such that
µk ≤ 1, for all k ≥ N1. (36) Fix any z ∈ Fix(T ). We claim
‖xk − z‖ ≤ max `∈[N1] {2R, ‖x` − z‖}, for all k ∈ N. (37)
The result (37) holds trivially for all k ∈ [N1]. Proceeding by induction, suppose (37) holds for some k > N1. If the inequality in Line 6 holds, then (35), (36) and the update formula in Line 7 together imply ‖xk+1‖ ≤ R. Since ‖S(z)‖ = 0, (35) also implies ‖z‖ ≤ R. Thus,
‖xk+1 − z‖ ≤ ‖xk+1‖+ ‖z‖ ≤ 2R ≤ max `∈[N1] {2R, ‖x` − z‖}. (38)
If instead the update in Line 9 is applied, the averagedness of T implies there is α ∈ (0, 1) such that
‖xk+1 − z‖2 ≤ ‖xk − z‖2 − 1− α α ‖S(xk)‖2 (39)
(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
‖xk+1 − z‖ ≤ ‖xk − z‖ ≤ max `∈[N1] {2R, ‖x` − z‖}. (40)
Equations (38) and (40) together close the induction, from which (37) follows. Whence
‖xk‖ ≤ ‖xk − z‖+ ‖z‖ ≤ max `∈[N1] {2R, ‖x` − z‖}+R, for all k ∈ N, (41)
which verifies the sequence {xk}k∈N is bounded.
Step 2: Because the inequality in Line 6 holds infinitely many times, there exists a subsequence {xqk}k∈N ⊆ {xk}k∈N satisfying
0 ≤ lim k→∞ ‖S(xqk)‖ ≤ lim k→∞ (1− δ)µk = 0, (42)
from which the squeeze theorem asserts ‖S(xqk)‖ → 0. Since {xk}k∈N is bounded, so also is {xqk}k∈N. Thus, there exists a subsequence {x`k}k∈N ⊆ {xqk}k∈N converging to a limit p ∈ H. Then applying the fact S is 2-Lipschitz and ‖ · ‖ is continuous yields
0 = lim k→∞ ‖S(x`k)‖ = ∥∥∥∥S ( limk→∞x`k )∥∥∥∥ = ‖S(p)‖ =⇒ p ∈ Fix(T ). (43) That is, {xk}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
Fix(T )[ε] := {x ∈ H : dFix(T )(x) ≤ ε}. (44)
Note Fix(T )[ε] is a nonempty closed and bounded subset ofH. Set C := ( 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. BecauseH is finite dimensional and C is 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 )[ε] =⇒ dFix(T )(x) ≤ ε, for all x ∈ H, (47)
where the first implication holds because x ∈ B(0, R) by (35) and x /∈ C by (46). By Assumption 3, there exists N2 ∈ N such that
µk ≤ ζ̃, for all k ≥ N2. (48) By the result of Step 2, there exists N3 ≥ N2 such that
‖xN3 − p‖ ≤ ε =⇒ dFix(T )(xN3) ≤ ε. (49) We claim
dFix(T )(x k) ≤ ε, for all k ≥ N3, (50)
which, by the arbitrariness of ε, implies (33) holds. Indeed, inductively suppose (50) holds for some k ≥ N3. If the inequality in Line 6 holds, then (47) and (48) together imply
‖S(xk+1)‖ ≤ (1− δ)µk ≤ ζ̃ =⇒ dFix(T )(xk+1) ≤ ε. (51)
Otherwise, letting P : H → H be the projection operator onto Fix(T ), we deduce
‖xk+1 − P (xk+1)‖ ≤ ‖xk+1 − P (xk)‖ ≤ ‖xk − P (xk)‖ = dFix(T )(xk) ≤ ε, (52)
where the second inequality follows from (39), taking z = P (xk). Note the left hand side of (52) equals the distance from xk+1 to Fix(T ). Therefore, (51) and (52) close the induction in each case, and (33) holds.
The limit (33) can be used in a similar manner to the work in Step 2 above to prove each cluster point of {xk}k∈N is in Fix(T ). Thus, if {xk}k∈N admits a unique cluster point, then the entire sequence converges to a point x? ∈ Fix(T ).
We restate Corollary 3.1 below and then provide a proof.
Corollary 3.1. If {xk}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.
Proof. The proof is parsed into four parts, one for each particular choice of the sequence {µk}k∈N in Table 2, where we note “Recent Term” is a special case of “Recent Max” obtained by taking m = 1. Each proof part is completely independent of the others and is separated by italic text. However, to avoid excessive writing, in each section let Γ ⊆ N be the set of all indices for which the descent condition in Line 6 holds, the sequence {tk}k∈N be an ascending enumeration of Γ, mk be the number of times the descent condition has been satisfied by iteration k, and µ1 ∈ (0,∞). Geometric Seqeunce. Define the sequence {µk}k∈N using, for each k ∈ N, the update formula
µk+1 =
{ µk if Line 6 holds,
(1− δ)µk otherwise. (53)
This implies µk = (1− δ)mkµ1. (54)
Since Γ is infinite, lim k→∞ mk =∞, and it follows that
lim k→∞ µk = lim k→∞
(1− δ)mkµ1 = 0 · µ1 = 0, (55)
i.e., Assumption 3 holds.
Arithmetic Average. Define the sequence {µk}k∈N using, for each k ∈ N, the update formula
µk+1 :=  1 mk + 1 ( ‖S(xk+1)‖+mkµk ) if Line 6 holds,
µk otherwise. (56)
Then observe
0 ≤ µtk+1 ≤ (1− δ)µtk +mtkµtk
mtk + 1 =
( 1− δ
mtk + 1
) µtk ≤ µtk , for all k ∈ N. (57)
Since µk+1 = µk whenever k /∈ Γ, (57) shows {µk}k∈N is monotonically decreasing. Consequently, using induction reveals
0 ≤ µtk − δ
mtk + 1 µtk ≤ µ1 − k∑ `=1 δµt` mt` + 1 = µ1 − k∑ `=1 δµt` `+ 1 for all k ∈ N, (58)
where we note mt` = ` in the sum since m` increments once each time a modification occurs in the sequence {µk}k∈N. By way of contradiction, suppose there exists α ∈ (0,∞) such that
lim inf k→∞
µk ≥ α > 0. (59)
Then (58) implies k∑ `=1 αδ `+ 1 ≤ k∑ `=1 δµt` `+ 1 ≤ µ1, for all k ∈ N. (60)
However, the sum on the left hand side becomes a divergent harmonic series as k →∞, contradicting the finite upper bound on the right hand side. This contradiction proves assumption (59) is false, from which it follows that
lim inf k→∞ µk = 0. (61)
By the monotone convergence theorem, we deduce µk → 0, i.e., Assumption 3 holds. Exponential Moving Average. Given θ ∈ (0, 1), for all k ∈ N, define
µk+1 :=
{ θ‖S(xk+1)‖+ (1− θ)µk−1 if Line 6 holds,
µk otherwise. (62)
Now observe
µtk+1 = θ‖S(xtk+1)‖+ (1− θ)µtk ≤ θ(1− δ)µtk + (1− θ)µtk = (1− θδ)µtk . (63)
This shows the sequence {µk}k∈N is nonincreasing and, when a decrease does occur, it is by a geometric factor of the current iterate. Through induction, it follows that
µk ≤ (1− θδ)mkµ1, for all k ∈ N. (64)
Since Γ is infinite, lim k→∞ mk =∞. This, combined with the fact (1− θδ) ∈ (0, 1), implies
0 ≤ lim k→∞ µk ≤ lim k→∞ (1− θδ)mkµ1 = 0 · µ1 = 0, (65)
from which the squeeze theorem asserts Assumption 3 holds.
Recent Max. Let m ∈ N. Set Ξk to be the set of the most recent min{m, k} indices for which the descent condition in Line 6 held, where {µk}k∈N is defined, for all k ∈ N, by the update formula
µk+1 =
{ max `∈Ξk ‖S(x`)‖ if Line 6 holds,
µk otherwise. (66)
The sequence {µk}k∈N is monotonically decreasing since the inequality in Line 6 implies, each time a new term ‖S(xk)‖ is introduced so that ‖S(xk)‖ ∈ Ξk+1, the new term is no larger than the largest term in Ξk. All that remains is to show this sequence converges to zero. By way of contradiction, suppose there exists α ∈ (0,∞) such that
lim inf k→∞ µk = α > 0. (67)
Then choose ε =
δα
2(1− δ) , (68)
which implies (1− δ)(α+ ε) < α. (69)
By (67) and the fact Γ is infinite, there exists Ñ ∈ N with Ñ > m such that
‖µtÑ − α| < ε =⇒ µtÑ < α+ ε. (70)
Then note each new element to Ξk is no larger than (1 − δ)µtÑ . And, for any k after m such replacements occur, it follows that
µk = max `∈Ξk
‖S(x`)‖ ≤ (1− δ)µtÑ ≤ (1− δ)(α+ ε) < α, (71)
a contradiction to (67). This contradiction shows our assumption (67) must be false, and so
lim inf k→∞ µk = 0. (72)
By the monotone convergence theorem, we conclude Assumption 3 holds.

Review summaries:
This paper proposes a framework to unfold the safeguarded Krasnosel’ski˘ı-Mann (SKM) method for the learn to optimization (L2O) schemes. First, SKM is proposed in Algorithm 1 with convergence guarantee established in Theorem 3.1 and Corollary 3.1. Then, SKM is unfolded and executed with a neural network summarized in Algorithm 2. Experiments on the Lasso and nonnegative least squares show the efficiency of the proposed method as well as the effectiveness of safeguarding compared to traditional L2O methods.   

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.