Distributionally-Robust Learning to Optimize Vinit Ranjan∗
Jisun Park†
Bartolomeo Stellato‡
arXiv:2605.06585v1 [cs.LG] 7 May 2026
May 8, 2026
Abstract We propose a distributionally robust approach to learning hyperparameters for first-order methods in convex optimization. Given a dataset of problem instances, we minimize a Wasserstein distributionally robust version of the performance estimation problem (PEP) over algorithm parameters such as step sizes. Our framework unifies two extremes: as the robustness radius vanishes, we recover classical learning to optimize (L2O); as it grows, we recover worst-case optimal algorithm design via PEP. We solve the resulting problem with stochastic gradient descent, differentiating through the solution of an inner semidefinite program at each step. We prove high-probability bounds showing that the true risk of the learned algorithm is at most the in-sample L2O optimum plus a slack that shrinks with the sample size, and is no worse than the worst-case PEP bound. On unconstrained quadratic minimization, LASSO, and linear programming benchmarks, our learned algorithms achieve strong out-of-sample performance with certifiable robustness, outperforming both worst-case optimal and vanilla L2O baselines.
1
Introduction
Applications in machine learning, signal processing, and optimal control increasingly rely on solving large-scale optimization problems using first-order methods for their low computational cost per iteration and modest memory requirements [9, 68]. Building on classical gradient descent [27], modern variants of first-order methods apply proximal operators [60] and operator splitting techniques [13, 15] to solve nonsmooth and constrained problems efficiently. This has led to robust open-source solvers such as PDLP [6, 7] for linear programs (LPs), OSQP [75] for quadratic programs (QPs), and SCS [59, 58], COSMO [28] for ∗
Department of Operations Research and Financial Engineering, Princeton University. Email: [email protected]. † Department of Operations Research and Financial Engineering, Princeton University, and Research Institute of Mathematics, Seoul National University. Email: [email protected]. ‡ Department of Operations Research and Financial Engineering, Princeton University. Email: [email protected].
1
semidefinite programs (SDPs) to name a few. Yet a fundamental challenge persists: practical performance depends critically on hyperparameters like step sizes, and tuning them well remains largely a manual effort prone to trial and error. Classical analysis derives worst-case convergence bounds over a given function class. The performance estimation problem (PEP) [24, 77] automates this approach by solving semidefinite programs to obtain tight worst-case guarantees. These bounds yield worst-case optimal constant step sizes [9], while recent work develops time-varying step sizes with improved worst-case rates [3]. However, worst-case analysis is inherently pessimistic: it guards against pathological instances that rarely arise in practice, yielding conservative algorithms that underperform on typical problems. Machine learning offers an alternative by exploiting structure in a dataset of problem instances. Learning to optimize (L2O) [5, 45, 18], also known as amortized optimization [4], unrolls algorithm iterations and tunes hyperparameters to minimize average loss on training data, often achieving substantial practical speedups over hand-tuned alternatives. However, L2O provides no theoretical guarantees on out-of-sample performance: learned algorithms may fail on problems outside the training distribution [18]. Recent efforts enforce guarantees via nonlinear systems theory [55, 56, 54] or PAC-Bayes bounds [25, 71, 76]. The former reintroduces worst-case conservatism through dynamical systems analysis; the latter yields bounds whose tightness depends sensitively on the prior and posterior used during training. In this work, we propose distributionally-robust learning to optimize (DR-L2O), a framework that combines computer-assisted worst-case analysis with data-driven L2O. Given a dataset of problem instances, we minimize the worst-case expected loss over algorithm parameters such as step sizes, where the worst case is taken over a Wasserstein ambiguity set centered at the empirical distribution. This builds on the recent distributionally robust PEP of [61], which evaluates this risk for a fixed algorithm; we instead minimize it, turning a performance certificate into a learning objective. A single Wasserstein radius controls the trade-off between data-driven and worst-case design: as the radius vanishes we recover L2O; as it grows we recover worst-case optimal design via PEP. We solve the resulting problem with stochastic gradient descent, where each iteration solves an inner SDP and differentiates through its solution via implicit differentiation of the KKT conditions. Figure 1 illustrates this trade-off. L2O minimizes the empirical loss but degrades on out-of-sample instances, the worst-case design is conservative on both, and DR-L2O sits between the two in-sample yet dominates both out-of-sample. Our contributions are as follows: 1. DR-L2O framework. We formulate the distributionally-robust learning to optimize problem (DR-L2O), which minimizes worst-case expected loss over a Wasserstein ambiguity set centered at the empirical distribution. We prove that varying the radius continuously interpolates between (L2O) and worst-case (OPT-PEP) (Proposition 2). 2. Scalable solution method. We solve (DR-L2O) via stochastic gradient descent, where each iteration solves an inner SDP and differentiates through its solution using implicit differentiation of the KKT conditions. 3. Out-of-sample guarantees. We prove that, for an appropriately chosen radius, 2
In-distribution
Out-of-distribution
Test loss
101
102
10−1
100
2
4
6
8
10
12
14
2
4
6
K L2O
8
10
12
14
K DR-L2O
OPT-PEP
Figure 1: Schematic comparison of L2O (L2O), the worst-case design (OPT-PEP), and (DR-L2O) on a fixed problem class. Solid lines show mean loss across instances; shaded bands show the 10th to 90th quantiles. Left: on the training distribution, L2O minimizes empirical loss, the worst-case design is conservative, and DR-L2O lies between (for most K). Right: on an out-of-distribution test set, L2O degrades with high variance, the worst-case design remains slow but robust, and DR-L2O dominates both for higher K. Details of this experiment are provided in Section 6.2.
the distributionally robust risk upper-bounds the true risk uniformly over algorithm parameters with high probability (Theorem 1). Combined with the interpolation result, the true risk of the learned algorithm is at most the in-sample L2O optimum plus a slack proportional to the radius, and at most the worst-case PEP bound (Theorem 2). 4. Numerical experiments. On unconstrained quadratic minimization, LASSO, and linear programs, we show that DR-L2O learns algorithms with strong out-of-sample performance and certifiable robustness, outperforming both worst-case optimal and vanilla L2O approaches.
1.1
Related works
Worst-case analysis and its limitations. The performance estimation problem (PEP) [24, 77] and control-theoretic analysis [44] automate worst-case convergence analysis for first-order methods in convex optimization, while verification frameworks [65, 64] extend these ideas to parametric quadratic and linear optimization. These approaches have enabled the discovery of worst-case optimal algorithms including the optimized gradient method (OGM) [39], OptISTA [36], the accelerated proximal point method (APPM) [38], and accelerated methods with silver step sizes [21, 3]. However, worst-case bounds are often far from typical performance. Average-case analysis [73, 80] offers tighter characterizations when problem distributions are known, while data-driven evaluation [72, 35, 61, 37] provides performance bounds on observed instances. Yet these approaches do not directly address how to learn hyperparameters that balance robustness with empirical performance. 3
Learning to optimize (L2O). L2O [5, 45, 18], also termed amortized optimization [4], learns algorithm hyperparameters from a dataset of problem instances. This has achieved practical success, from learned proximal methods that solve LASSO in only a few iterations [47] to scalable approaches for neural network training [19]. However, L2O is essentially empirical risk minimization, offering no guarantees on out-of-sample performance. Two strategies have emerged to address this gap. The first imposes structural constraints that ensure convergence: [20] identified algorithm structures that are both convergent and trainable, [63] and [33] developed safeguarded formulations guaranteeing convergence for every instance, and [48, 55, 56] derived necessary conditions for convergent learned algorithms. While effective, these structural approaches limit learning capacity. The second strategy incorporates distribution-aware objectives: [81] proposed meta-learning for better generalization, [74] developed out-of-distribution-robust training via data alignment, [71, 76] employed PAC-Bayes bounds, and [70] used worst-case convergence rates as a regularizer in the learning objective. Both strategies have limits: structural constraints reduce learning capacity, while distribution-aware objectives either reintroduce worst-case conservatism or rely sensitively on PAC-Bayes priors and posteriors. Distributionally robust optimization (DRO). DRO [42] handles uncertainty by optimizing over ambiguity sets, which are collections of distributions consistent with partial information. When only data samples are available, data-driven DRO [79, 42] constructs ambiguity sets as metric balls centered at the empirical distribution, using Wasserstein distance [57, 41] or divergences such as Kullback-Leibler (KL) divergence [34]. This framework naturally quantifies distributional uncertainty from finite samples, providing the mechanism we exploit to interpolate between L2O and worst-case PEP.
2
Problem setup
We consider a problem instance z = (f, x0 ) consisting of a convex optimization problem minimize f (x), x
where f : Rd → R ∪ {∞} is convex, proper, and lower semi-continuous, and x0 ∈ Rd is an initial point. We assume that the problem has a solution x⋆ with optimal value f ⋆ = f (x⋆ ). An algorithm Aθ maps a problem instance to a sequence of approximate solutions Aθ (z) = {xk }k=0,1,... , where θ ∈ Θ represents the algorithm parameters such as step sizes. We restrict attention to first-order algorithms, which use only subgradient information of f . Assumption 1 (First-order algorithms). For any z = (f, x0 ), the sequence {xk }k=0,1,... = Aθ (z) satisfies the span condition xk ∈ x0 + span{g 0 , . . . , g k−1 , g k },
4
k = 0, 1, . . . ,
where g i ∈ ∂f (xi ) for i = 0, . . . , k. Note xk is ∂f (xk )-dependent when the algorithm involves proximal step xk = proxηf (xk−1 ), i.e., g k = xk−1 − xk ∈ ∂f (xk ). We measure algorithm performance on instance z by a loss function ℓ(Aθ (z)) ∈ R+ . Given an iteration budget K, common choices of ℓ are the function-value gap f (xK ) − f ⋆ and the distance to optimality ∥xK − x⋆ ∥2 , evaluated at the K-th iterate xK . Let P be a probability distribution over the set Z of problem instances. Our goal is to find the optimal algorithm parameter θ⋆ minimizing the risk, defined as the expected loss under P: minimize R(θ, P) := E ℓ(Aθ (z)) . (1) z∼P θ∈Θ
3
Performance guarantees
We now describe two approaches to measuring algorithm performance that serve as building blocks for our framework: worst-case guarantees via PEP [24, 77] and data-driven probabilistic guarantees via Wasserstein DRO [61].
3.1
Worst-case guarantee
The worst-case loss supz∈Z ℓ(Aθ (z)) upper-bounds the risk R(θ, P) for any distribution P supported on Z. It is evaluated as the optimal value of the following performance estimation problem (PEP): maximize ℓ(Aθ (z)) (PEP) subject to z = (f, x0 ) ∈ Z = F × X ,
where F is the function class and X is the set of initial iterates. When the function class F admits interpolation conditions, e.g., smooth convex or convex Lipschitz functions, (PEP) admits a tractable formulation as a convex conic program [24, 77, 67]. Following [24, 77], (PEP) admits an equivalent SDP via a lifted representation. Given K+2 0 T the iterates {xk }K with k=0 = Aθ (f, x ), we define the Gram matrix G = P P ∈ S+ h i P = x0 − x⋆ g 0 · · · g K ∈ Rd×(K+2) where g k ∈ ∂f (xk ), and the function-value vector F = (f (x0 ) − f ⋆ , . . . , f (xK ) − f ⋆ ) ∈ RK+1 . Each component of (PEP) is linear in (G, F ): the loss is ℓ(G, F ) = tr(ATobj G) + bTobj F , the initial condition is tr(AT0 G) + bT0 F + c0 ≤ 0, and the interpolation constraints of F are Sθ (G, F ) ∈ RM + with Sθ (G, F ) m = − tr Am (θ)T G − bm (θ)T F, m = 1, . . . , M, (2)
where the LMI coefficients Am (θ) ∈ SK+2 and bm (θ) ∈ RK+1 depend on the algorithm parameters θ [37]. We denote by Ξθ the feasible set of (G, F ) pairs: T T Ξθ = (G, F ) ∈ SK+2 × RK+1 Sθ (G, F ) ∈ RM + + , tr(A0 G) + b0 F + c0 ≤ 0 . The explicit SDP and its Lagrangian dual are given in Section A.1. 5
3.2
Data-driven probabilistic guarantee
Since the true distribution P is unknown, we approximate the risk in (1) from a finite dataset. PN N b bN = {ẑi }N Given i.i.d. samples D i=1 ∼ P , the empirical distribution is PN = (1/N ) i=1 δẑi , where δz is the Dirac measure at z ∈ Z. The corresponding empirical risk is b N ) = (1/N )PN ℓ Aθ (ẑi ) . R(θ, P i=1 With finite samples, however, the empirical risk may be a poor estimate of the true bN yields a lifted representation risk R(θ, P). Running Aθ on each sample ẑi ∈ D θ b θ = (1/N ) PN δ b b over Ξθ . bi , Fbi ) ∈ Ξ , inducing the lifted empirical distribution P (G N i=1 (Gi ,Fi ) To mitigate the bias of the empirical risk, [61] proposes a data-driven framework that yields a probabilistic performance guarantee robust to sampling variability. The guarantees are derived from a Wasserstein distributionally robust variant of the performance estimation problem (DRO-PEP): b N ) := maximize R(θ, Q) Rε (θ, P b θ ), subject to Q ∈ Uε (P N
(DRO-PEP)
b θ ) is the 1-Wasserstein ball of radius ε centered at P bθ : where the ambiguity set Uε (P N N n o b θ ) = Qθ supp Qθ ⊆ Ξθ , W1 (P b θ , Qθ ) ≤ ε , Uε (P N N with W1 induced by the norm ∥(G, F )∥ = DRO risk.
p b N ) as the ∥G∥2F + ∥F ∥2 . We refer to Rε (θ, P
Tractable formulation. By [61, Theorem 2], (DRO-PEP) has a tractable formulation as a convex conic program; see Section A.2 for details. For θ ∈ Θ and ε > 0, (DRO-PEP) also admits the following equivalent robust optimization formulation: P T T maximize (1/N ) N i=1 tr(Aobj Gi ) + bobj Fi PN bi , Fbi )∥ ≤ ε subject to (1/N ) i=1 ∥(Gi , Fi ) − (G (Gi , Fi ) ∈ Ξθ , i = 1, . . . , N.
(DRO-PEP-P)
The equivalence comes from the strong duality between the dual of (DRO-PEP) and its bi-dual, which is exactly (DRO-PEP-P). The proof of the strong duality result is deferred to Section A.3. Proposition 1. Let θ ∈ Θ and ε > 0. Then (DRO-PEP-P) is a bi-dual of (DRO-PEP), and their optimal values coincide. Furthermore, (DRO-PEP-P) has an optimal solution. The key departure of this work from [61] is to minimize, rather than evaluate, the DRO risk over θ, turning a performance certificate into a learning objective. 6
4
Distributionally-robust learning to optimize
b N ) over θ to learn optimal algorithm hyperparameWe now minimize the DRO risk Rε (θ, P ters. Unlike [61], which uses the DRO risk only to certify performance of a fixed θ, we treat it as a learning objective. When the learning objective is the worst-case loss, the corresponding learning problem is [37] minimize sup ℓ Aθ (z) . (OPT-PEP) θ∈Θ
z∈Z
However, worst-case optimal algorithms can underperform on typical instances. At the other extreme, learning to optimize (L2O) minimizes the empirical risk [18]: b N ) = (1/N ) minimize R(θ, P θ∈Θ
PN
ℓ A (ẑ ) . θ i i=1
(L2O)
When the loss ℓ is differentiable, (L2O) is solved efficiently with stochastic first-order methods. However, it provides no out-of-sample generalization guarantee: the learned algorithm Aθ⋆ may fail on unseen instances. To balance the two extremes, we minimize the DRO risk, yielding the learning problem b N ), minimize Rε (θ, P θ∈Θ
(DR-L2O)
which we call distributionally-robust learning to optimize.
4.1
Solution method
bi , Fbi , Am (θ), and bm (θ) are nonProblem (DR-L2O) is non-convex in θ, since the entries of G b N ) is the optimal value of convex quadratic polynomials in θ. Still, the DRO risk Rε (θ, P b b N the convex conic program (DRO-PEP-D) parametrized by {(Am , bm )}M m=1 and {(Gi , Fi )}i=1 , whose solution map is differentiable using the approach of [2]. We compute the gradib N )/dθ via the chain rule and apply stochastic gradient descent to obtain θ⋆ . ent dRε (θ, P See Figure 2 for a summary of the solution method and Section C for a detailed description of the gradient evaluation and solution algorithm; the latter also establishes that (DR-L2O) attains its minimum over Θ (Remark 3).
5
Theoretical guarantees of learned optimizer
In this section, we show that (DR-L2O) not only interpolates between (OPT-PEP) and (L2O), but also yields out-of-sample and out-of-distribution guarantees. Throughout this section, we assume that the parameter set Θ is compact. This holds in practice because of finite floating-point precision. Assumption 2. The parameter set Θ is compact, i.e., closed and bounded.
7
! <latexit sha1_base64="nF4MVGkFocKxmGy8t9/fobrr/CI=">AAACC3icbZDLSsNAFIYnXmu9RV26GVqEClISkeqy6MZlBXuBJpbJdNIOnVyYORFL6N6Nr+LGhSJufQF3vo2TNgttPTDw8f/ncOb8Xiy4Asv6NpaWV1bX1gsbxc2t7Z1dc2+/paJEUtakkYhkxyOKCR6yJnAQrBNLRgJPsLY3usr89j2TikfhLYxj5gZkEHKfUwJa6pklx+ODCnaGBFJ/0uMnM3zQeGfhzDzumWWrak0LL4KdQxnl1eiZX04/oknAQqCCKNW1rRjclEjgVLBJ0UkUiwkdkQHragxJwJSbTm+Z4COt9LEfSf1CwFP190RKAqXGgac7AwJDNe9l4n9eNwH/wk15GCfAQjpb5CcCQ4SzYHCfS0ZBjDUQKrn+K6ZDIgkFHV9Rh2DPn7wIrdOqXavWbs7K9cs8jgI6RCVUQTY6R3V0jRqoiSh6RM/oFb0ZT8aL8W58zFqXjHzmAP0p4/MHkViaIQ==</latexit>
fˆi , x̂0i
"
problem instance
Algorithm
Aθ
<latexit sha1_base64="TJzjEHIgMGCR2jZ9D0ra+G+jFuU=">AAACFXicdVBNa9tAEF2lzUedjzrtsZelJpBDECsTu/EttJce0xInAcuY0WoUL1mtxO4oYIT+RC75K73kkFJ6LfTWf9O140BT2gfLPt6bYWZeUmrlSIhfwcqz56tr6xsvWptb2zsv27uvzlxRWYlDWejCXiTgUCuDQ1Kk8aK0CHmi8Ty5+jD3z6/ROlWYU5qVOM7h0qhMSSAvTdoHcWZB1imPc6CpBF1/biZ1fA0WS6d0YZqmTmOaIkHDJ+2OCEX3qCcGXIRd/3V7nvRENOgPeBSKBTpsiZNJ+2ecFrLK0ZDU4NwoEiWNa7CkpMamFVcOS5BXcIkjTw3k6Mb14qqG73kl5Vlh/TPEF+qfHTXkzs3yxFfOd3d/e3PxX96oouxoXCtTVoRGPgzKKs2p4POIeKosStIzT0Ba5Xflcgo+JvJBtnwIj5fy/5Ozbhj1w/6nw87x+2UcG+wNe8v2WcTesWP2kZ2wIZPshn1h9+xrcBvcBd+C7w+lK8Gy5zV7guDHb+C5oJM=</latexit>
dRε dθ
b i , Fbi ) (G <latexit sha1_base64="dl5EEuXQ/UNVfBb2KmjCAjIrT/k=">AAACB3icbVDLSgMxFM3UV62vUZeCBItQQcqMSHVZFNRlBfuAdhgymUwbmnmQ3FHK0J0bf8WNC0Xc+gvu/BvTB6LVA4Fzz7mXm3u8RHAFlvVp5ObmFxaX8suFldW19Q1zc6uh4lRSVqexiGXLI4oJHrE6cBCslUhGQk+wptc/H/nNWyYVj6MbGCTMCUk34gGnBLTkmrulzh33WY9Adjl0+SH+Li90eeCaRatsjYH/EntKimiKmmt+dPyYpiGLgAqiVNu2EnAyIoFTwYaFTqpYQmifdFlb04iETDnZ+I4h3teKj4NY6hcBHqs/JzISKjUIPd0ZEuipWW8k/ue1UwhOnYxHSQosopNFQSowxHgUCva5ZBTEQBNCJdd/xbRHJKGgoyvoEOzZk/+SxlHZrpQr18fF6tk0jjzaQXuohGx0gqroCtVQHVF0jx7RM3oxHown49V4m7TmjOnMNvoF4/0Lc4KZEA==</latexit>
<latexit sha1_base64="0bSqywf1v009Ih6g0tEnMMFFOPo=">AAAB/XicbVDLSsNAFJ34rPUVHzs3g0VwVRKR6rLqxmUF+4AmhMl00g6dPJi5EWoI/oobF4q49T/c+TdO2iy09cDA4Zx7uWeOnwiuwLK+jaXlldW19cpGdXNre2fX3NvvqDiVlLVpLGLZ84ligkesDRwE6yWSkdAXrOuPbwq/+8Ck4nF0D5OEuSEZRjzglICWPPPQCQmMKBHZVe5lDowYkNwza1bdmgIvErskNVSi5ZlfziCmacgioIIo1betBNyMSOBUsLzqpIolhI7JkPU1jUjIlJtN0+f4RCsDHMRSvwjwVP29kZFQqUno68kiq5r3CvE/r59CcOlmPEpSYBGdHQpSgSHGRRV4wCWjICaaECq5zorpiEhCQRdW1SXY819eJJ2zut2oN+7Oa83rso4KOkLH6BTZ6AI10S1qoTai6BE9o1f0ZjwZL8a78TEbXTLKnQP0B8bnDy0MlbY=</latexit>
simulate trajectories
bN) Rε (θ, P <latexit sha1_base64="MBf6IsUN54tr5UxyTAN5djCwf1I=">AAACI3icbVBBSxtBGJ1Va21aNdVjL4NBSEHCrkhaepL20lNJxaiQDcu3k2/NkNnZZeZbS1j2v/TiX/HSQ0V66aH/pbNrDhp9MPB4733M9704V9KS7//1VlbXXqy/3HjVev1mc2u7/XbnzGaFETgUmcrMRQwWldQ4JEkKL3KDkMYKz+PZl9o/v0JjZaZPaZ7jOIVLLRMpgJwUtT+FKdBUgCpPqqgMr8BgbqXKdMW7IU2R4ICHP+QEp0Blk42TclBV0bf3Ubvj9/wG/CkJFqTDFhhE7btwkokiRU1CgbWjwM9pXIIhKRRWrbCwmIOYwSWOHNWQoh2XzY0V33fKhCeZcU8Tb9SHEyWk1s7T2CXrLe2yV4vPeaOCko/jUuq8INTi/qOkUJwyXhfGJ9KgIDV3BISRblcupmBAkKu15UoIlk9+Ss4Oe0G/1/9+1Dn+vKhjg71je6zLAvaBHbOvbMCGTLCf7Ib9ZrfetffLu/P+3EdXvMXMLnsE799/002lnA==</latexit>
distributionally-robust performance estimation
DRO risk
<latexit sha1_base64="UxUKi+81CzMKEo2qh7SU2nR5SyY=">AAACJHicbVDLSgMxFM34rOOr6tJNsAi6KTMiVXBTdeOygtVCp5RMetsGM5khuSOUoR/jxl9x48IHLtz4LabtgFo9EDicc26Se8JECoOe9+HMzM7NLywWltzlldW19eLG5rWJU82hzmMZ60bIDEihoI4CJTQSDSwKJdyEt+cj/+YOtBGxusJBAq2I9ZToCs7QSu3iSRBCT6iMac0Gw4wPXUpP29FegH1Atk+DwArht+AGoDp5ul0seWVvDPqX+DkpkRy1dvE16MQ8jUAhl8yYpu8l2LK3oeAShm6QGkgYv2U9aFqqWASmlY2XHNJdq3RoN9b2KKRj9edExiJjBlFokxHDvpn2RuJ/XjPF7nErEypJERSfPNRNJcWYjhqjHaGBoxxYwrgW9q+U95lmHG2vri3Bn175L7k+KPuVcuXysFQ9y+sokG2yQ/aIT45IlVyQGqkTTu7JI3kmL86D8+S8Oe+T6IyTz2yRX3A+vwDYeKO4</latexit>
Am (θ) bm (θ)
interpolation conditions
<latexit sha1_base64="38yPlnNU7XSqoFQtZv9wi5VhQAI=">AAACO3icdVDPSxtBFJ61trXpr7QeexkMgoWwzAaTGnrRFqpHlUaFbFhmZ9+awZndZeatEpb9v3rpP9FbL7300CJevTuJqWjRB8N873vfx8z74kJJi4z99BYeLT5+8nTpWeP5i5evXjffvD2weWkEDESucnMUcwtKZjBAiQqOCgNcxwoO45PP0/nhKRgr8+wrTgoYaX6cyVQKjo6KmvtharioEroWnskExhyr7TqSbXrTfnHt+9opQhwD8rodfmzcmLYi3aZxpG8LomaL+ayz0WV9yvyOuzpdB7os6Pf6NPDZrFpkXrtR80eY5KLUkKFQ3NphwAocVdygFArqRlhaKLg44ccwdDDjGuyomu1e01XHJDTNjTsZ0hl721Fxbe1Ex06pOY7t/7Mped9sWGK6MapkVpQImbh+KC0VxZxOg6SJNCBQTRzgwkj3VyrG3OWCLu6GC+HfpvRhcNDxg57f21tvbX6ax7FE3pEVskYC8oFskh2ySwZEkG/kF/lD/nrfvd/euXdxLV3w5p5lcqe8yyuEJK3O</latexit>
b i , Fbi ) d(Am , bm ) d(G , dθ dθ
<latexit sha1_base64="azNUSXYbsEKQwGP5SEqN6OfWrl8=">AAACfnicnVFLj9MwEHbCaymvLhy5GCrQLpTgVNvSissCEnBcEN1dqamiiTvZWms7ke0sqqL8DP4YN34LF5xueQoujGT5m29mPJ5vslIK6xj7EoQXLl66fGXraufa9Rs3b3W3bx/aojIcp7yQhTnOwKIUGqdOOInHpUFQmcSj7PRVGz86Q2NFoT+4VYlzBSda5IKD81Ta/ZTkBnidlGCcAEkTBW7JQdbvm7ROzsBgaYUsdNP8zNl5kao+zVK12/ST553/eSH5KBa4BFe/aVLRpz/c197dbTppt8ciNhgP2YSyaOCvwdCDIYsnowmNI7a2HtnYQdr9nCwKXinUjkuwdhaz0s3rth2X2HSSymIJ/BROcOahBoV2Xq/la+gDzyxoXhh/tKNr9teKGpS1K5X5zHY2+2esJf8Wm1UuH89rocvKoebnjfJKUlfQdhd0IQxyJ1ceADfC/5XyJXg1nd9YK8L3Sem/weEgikfR6N1eb//lRo4tcpfcJzskJs/IPnlLDsiUcPI1uBc8Ch6HJHwYPgmfnqeGwabmDvnNwvE3Qj7EGg==</latexit>
∂Rε ∂Rε , b i , Fbi ) ∂(Am , bm ) ∂(G
Figure 2: The diagram of solution method solving (DR-L2O) with stochastic gradient methods. Given θ, we have interpolation LMI coefficients Am (θ) and bm (θ), while a sample problem inb i , Fbi ) = Aθ (fˆi , x̂0 ). stance (fˆi , x̂0i ) is mapped to the grammian and function-value representation (G i b N ) w.r.t. θ through back-propagation using chain We evaluate the gradient of DRO risk Rε (θ, P rule.
5.1
Robustness certificate
We first show that our learned algorithm enjoys a robust out-of-sample guarantee, a consequence of the data-driven Wasserstein DRO framework [57]. Here, we extend the finitesample guarantee of [61, Section 3.2] from a fixed θ ∈ Θ to a uniform statement over θ ∈ Θ. The proof is provided in Section B.1. Theorem 1. Suppose Assumption 2 holds. Let N and K be fixed. For each β ∈ (0, 1), there exists ε(β) > 0 and L > 0 such that ! PN E ℓ(Aθ (z)) ≤ E sup ℓ(G, F ) , ∀θ ∈ Θ ≥ 1 − β, z∼P
bθ ) Qθ ∈Uε(β) (P N
(G,F )∼Qθ
where q is the dimension of (G, F ) and ε(β) ≲
log(1/β) + dim(Θ) log N + log diam(Θ) N
!1/q +
2L . N
Remark 1 (Out-of-distribution (OOD) guarantee). The DRO risk provides a performance guarantee that is also robust to the choice of true distribution, in the sense that the same b θ⋆ ), not just the true distribution. bound holds for any Qθ ∈ Uε (P N
5.2
Interpolating behavior of (DR-L2O)
We now show how (DR-L2O) interpolates between (L2O) and the worst-case optimal design (OPT-PEP) as the radius ε varies. All proofs are presented in Section B.2. 8
Proposition 2. Suppose Assumption 2 holds, ε > 0, and let θε⋆ be an optimal solution b N ) is monotonically nondecreasing. Furtherof (DR-L2O). Then the mapping ε 7→ Rε (θε⋆ , P more, b N ) = inf R(θ, P b N ), b N ) = inf sup ℓ Aθ (z) . lim+ Rε (θε⋆ , P lim Rε (θε⋆ , P ε→∞
θ∈Θ
ε→0
θ∈Θ z∈Z
As ε → 0+ , problem (DR-L2O) loses its distributional robustness and reduces to (L2O). For large enough ε, the optimizer learned via (DR-L2O) is robust to essentially all distributions supported on Ξ, becoming worst-case optimal. We now establish certified performance bounds for the true risk of our learned algorithm. Theorem 2. Suppose Assumption 2 holds. Given ε > 0, let θε⋆ be an optimal solution of (DR-L2O). For any β ∈ (0, 1) and ε = ε(β) as in Theorem 1, the true risk R(θε⋆ , P) of the optimizer learned via (DR-L2O) satisfies the following with probability at least 1 − β: b N ) + ε Lip(ℓ), R(θε⋆ , P) ≤ inf R(θ, P θ∈Θ
R(θε⋆ , P) ≤ inf sup ℓ Aθ (z) , θ∈Θ z∈Z
where Lip(ℓ) is a Lipschitz constant of the loss function ℓ. The first bound certifies that the true risk of the optimizer learned via (DR-L2O) is within ε Lip(ℓ) of the in-sample optimum of (L2O), providing an out-of-sample guarantee that (L2O) alone lacks. The second bound validates that the true risk is no worse than the worst-case optimal bound from (OPT-PEP). Note that these bounds do not imply strict empirical superiority over either baseline; the numerical experiments in Section 6 demonstrate that in practice (DR-L2O) achieves lower true risk than both.
6
Numerical experiments
We now showcase the performance of our learned optimizerson unconstrained quadratic minimization, LASSO, and total variation inpainting. For each experiment, we report the out-of-sample test loss, the fraction of test instances solved at relative tolerance η > 0 (an instance z = (f, x0 ) is solved if ℓ(Aθ (z)) ≤ η(1 + |f ⋆ |)), and out-of-distribution performance. Implementation details and per-iteration wall-time tables are in Appendix D. We compare ⋆ three learned optimizers with parameters θε⋆ from (DR-L2O), θPEP from (OPT-PEP), and ⋆ θL2O from (L2O). All code to reproduce our experiments is available at https://github.com/stellatogrp/dro_pep. Training and parameter selection. We learn a step-size schedule θ = {θk }K−1 k=0 for each horizon K and each of (DR-L2O), (OPT-PEP), and (L2O) using gradient-based methods. We cross-validate the Wasserstein radius ε, the learning rate, and the AdamW weight decay, and use a held-out test set for final loss evaluations. Full details are in Section C. 9
Test set, fraction of problems solved
Out-of-distribution
In-distribution
η = 0.001
η = 0.01
η = 0.1
1
1
1
0
0
0
1
1
1
0
0 2
4
6
8 10 12 14 K
0 2
L2O
4
6
8 10 12 14 K
DR-L2O
2
4
6
8 10 12 14 K
OPT-PEP
Figure 3: Quadratic minimization experiment fractions of problems solved to different relative tolerances across horizon K. Top: Fraction solved on a held-out in-distribution test set. Bottom: Fraction solved on an out-of-distribution set.
6.1
Unconstrained quadratic minimization
Consider an unconstrained quadratic minimization problem minimize (1/2)xT Qx, x
K where Q ∈ Sd++ is a positive definite matrix. We learn the step size θ = (θk )K−1 k=0 ∈ R+ of gradient descent xk+1 = (I − θk Q)xk for k = 0, . . . , K − 1. For the performance loss, we use the objective value gap f (xK ) − f (x⋆ ), where f (x⋆ ) = 0 for this quadratic function class. We sample Qi from the Marčenko-Pastur distribution [53, 62] and zi0 uniformly on a ball; the out-of-distribution shift increases the Lipschitz constant L (full parameters in Appendix D.1).
Results. In Figure 3 we show the fractions of problems solved by the learned schedules. For the in-distribution test set, (L2O) solves the most problems at the tolerance levels but (DR-L2O) performs the best out-of-distribution for larger K. Figure 7 (Appendix D.1) shows that the (DR-L2O) loss remains competitive for many instances but the average degrades with K due to some instances where the learned schedule performs poorly. Both (OPT-PEP) and our (DR-L2O) are robust against the distribution shift while the (DR-L2O) average performance does not degrade.
10
Test set, fraction of problems solved
Out-of-distribution
In-distribution
η = 0.01
η = 0.05
η = 0.1
1
1
1
0
0
0
1
1
1
0
0 2
4
6
8 10 12 14 K
0 2
L2O-ISTA
4
6
8 10 12 14 K
L2O-ALISTA
DR-L2O
2
4
6
8 10 12 14 K
OPT-PEP
Figure 4: Fractions of LASSO problems solved to different relative tolerances across horizon K. Top: Fraction solved in-distribution. Bottom: Fraction solved out-of-distribution.
6.2
LASSO
Consider an ℓ1 -regularized least-squares problem: minimize (1/2)∥Ax − b∥2 + λ∥x∥1 , x
where λ > 0 is the ℓ1 regularization parameter. We learn the step sizes {θk }K−1 k=0 of ISTA [9]: xk+1 = proxλθk ∥·∥1 xk − θk AT (Axk − b) , K ⋆ K where θ = (θk )K−1 k=0 ∈ R+ . We again consider the performance loss f (x ) − f (x ). Specialized learned optimizers exist for LASSO, all based on empirical risk minimization. Learned ISTA (LISTA) [31] learns weight matrices to speed up the proximal updates over a distribution of instances. For comparison, we label our default (L2O) scheme as L2O-ISTA and compare against analytic LISTA (ALISTA) [20, 47], which has better practical performance and is more computationally efficient than LISTA. ALISTA restricts the structure of LISTA by making the update rule xk+1 = proxγ k ∥·∥1 xk − θk W T (Axk − b) , where W is precomputed as the solution of a convex quadratic program [20, Equation 16]. We label this baseline L2O-ALISTA. We use a sparse coding setup, recovering sparse vectors x̃ from noisy measurements b = Ax̃ + ε with a fixed dictionary A ∈ Rm×n [20]; the out-of-distribution shift increases the sparse-signal variance σx (full parameters in Appendix D.2).
Results. We previewed the loss results of this experiment in Figure 1, and the full version with L2O-ALISTA is provided in Figure 8 (Appendix D.2), while the fractions of problems solved are provided in Figure 4. Overall, we notice similar trends as in the unconstrained 11
quadratic experiment. The biggest difference is that, even for the same horizon K, the performance of L2O-ISTA degrades significantly more. Given that the 10th quantile is so low and the schedule is able to solve many of the out-of-distribution problems, this implies that the schedule diverges on a small number of instances, while our (DR-L2O) learned schedule is robust against such outliers.
6.3
Total variation (TV) inpainting
We consider the ℓ1 total-variation (TV) inpainting problem, which fills in an image that has been corrupted with some pixels blacked out [23]. For a grayscale image represented as an m × n matrix U orig with values in [0, 255] (scaled down to [0, 1] for algorithm stability) and the set K of known pixel locations, the goal is to reconstruct the image U while matching the known pixels. The corresponding linear program is " # m−1 n−1 XX Ui+1,j − Ui,j minimize U Ui,j+1 − Ui,j i=1 j=1 1
subject to Uij = Uijorig , Uij ∈ [0, 1],
(i, j) ∈ K i = 1, . . . , m, j = 1, . . . , n.
We solve this LP using primal-dual hybrid gradient (PDHG) [15] adjusted for linear programs (PDLP) [7]. We first cast the problem in standard LP form by introducing slack variables for the ℓ1 terms (see Section D.3), and apply PDHG to its convex-concave saddle reformulation minimize maximize L(x, u) = f (x) + ⟨u, M x⟩ − g ∗ (u), x
u
where f encodes the LP’s linear cost and box constraints, M stacks the equality and inequality constraint matrices, and g(M x) is the indicator of the constraints; see Section D.3 for the explicit form and the corresponding PDHG iteration, whose step sizes 3K θ = {(τ k , ρk , σ k )}K−1 k=0 ∈ R+ we learn. For the performance loss, we use the Lagrangian duality gap L(x, u⋆ ) − L(x⋆ , u); since this class falls outside the standard (PEP) interpolation conditions, we extend them to general linear operators in Section A.4. For in-distribution we use Olivetti faces [69, 14] with 10% pixels blacked out; for out-of-distribution we use color images from Tiny ImageNet [43, 22] (full parameters in Appendix D.3). Results. We show both the fractions of problems solved and the losses in Figures 10 and 9 in Appendix D.3. Here, in Figure 5, we instead show the reconstructions of images from the two datasets with the different learned schedules. The Olivetti reconstruction is not significantly different between (L2O) and (DR-L2O). However, the Tiny ImageNet reconstruction shows a difference where the learned schedules increasingly bring the image into focus and lessen the gray hue. In Figure 11 (Appendix D.3), we provide more reconstructions on extra images from Tiny ImageNet.
12
Corrupted
Opt. LP Reconstruction
OPT-PEP
L2O
DR-L2O
Tiny-ImageNet
Olivetti
Original
Figure 5: TV inpainting reconstructions (K = 10). Left to right: original, corrupted (10% missing pixels), optimal LP reconstruction, and reconstructions from (OPT-PEP), (L2O), and (DR-L2O).
7
Conclusion
We presented a distributionally-robust algorithm design framework for first-order methods. By minimizing the distributionally robust risk, our learned algorithms enjoy out-of-sample and out-of-distribution guarantees that classical L2O approaches lack. The framework extends naturally to risk-averse training with other risk measures, such as conditional valueat-risk [66]; the corresponding DRO-PEP formulation already appears in [61, Corollary 2.1]. A major limitation is scalability: unrolling algorithms for more iterations requires solving a much larger conic program at each gradient step. A natural remedy is to use GPU-accelerated first-order solvers [52, 46, 51, 16, 17].
13
Acknowledgements Bartolomeo Stellato and Vinit Ranjan are supported by the NSF CAREER Award ECCS2239771 and the ONR YIP Award N000142512147. Jisun Park is supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (RS2024-00353014) and the ONR YIP Award N000142512147. The authors are pleased to acknowledge that the work reported on in this paper was substantially performed using Princeton University’s Research Computing resources.
References [1] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter. Differentiable Convex Optimization Layers. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. [2] A. Agrawal, S. Barratt, S. Boyd, E. Busseti, and W. M. Moursi. Differentiating Through a Cone Program. Journal of Applied & Numerical Optimization, 1(2):107–115, May 2019. [3] J. M. Altschuler and P. A. Parrilo. Acceleration by stepsize hedging: Silver Stepsize Schedule for smooth convex optimization. Mathematical Programming, 213(1):1105– 1118, Sept. 2025. [4] B. Amos. Tutorial on Amortized Optimization. Foundations and Trends in Machine Learning, 16(5):592–732, June 2023. [5] M. Andrychowicz, M. Denil, S. Gómez, M. W. Hoffman, D. Pfau, T. Schaul, B. Shillingford, and N. de Freitas. Learning to learn by gradient descent by gradient descent. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016. [6] D. Applegate, M. Diaz, O. Hinder, H. Lu, M. Lubin, B. O’ Donoghue, and W. Schudy. Practical Large-Scale Linear Programming using Primal-Dual Hybrid Gradient. In Advances in Neural Information Processing Systems, volume 34, pages 20243–20257, New York, 2021. Curran Associates, Inc. [7] D. Applegate, O. Hinder, H. Lu, and M. Lubin. Faster first-order primal-dual methods for linear programming using restarts and sharpness. Mathematical Programming, 201(1):133–184, 2023. [8] H. H. Bauschke and P. L. Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. CMS Books in Mathematics. Springer International Publishing, Cham, 2017.
14
[9] A. Beck. First-Order Methods in Optimization. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017. [10] C. Berge. Topological Spaces: Including a Treatment of Multi-valued Functions, Vector Spaces and Convexity. Oliver & Boyd, 1963. [11] N. Blin, S. Gualandi, C. Maes, A. Lodi, and B. Stellato. Batched First-Order Methods for Parallel LP Solving in MIP. In International Conference on Machine Learning (ICML), 2026. [12] N. Bousselmi, J. M. Hendrickx, and F. Glineur. Interpolation Conditions for Linear Operators and Applications to Performance Estimation Problems. SIAM Journal on Optimization, 34(3):3033–3063, Sept. 2024. [13] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning, 3:1–122, Jan. 2011. [14] A. L. Cambridge. The Olivetti faces dataset. http://www.cl.cam.ac.uk/research/ dtg/attarchive/facedatabase.html, 1994. [15] A. Chambolle and T. Pock. A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging. Journal of Mathematical Imaging and Vision, 40(1):120– 145, 2011. [16] K. Chen, D. Sun, Y. Yuan, G. Zhang, and X. Zhao. HPR-LP: An implementation of an HPR method for solving linear programming. Mathematical Programming Computation, Oct. 2025. [17] K. Chen, D. Sun, Y. Yuan, G. Zhang, and X. Zhao. HPR-QP: A dual Halpern PeacemanRachford method for solving large-scale convex composite quadratic programming, July 2025. [18] T. Chen, X. Chen, W. Chen, H. Heaton, J. Liu, Z. Wang, and W. Yin. Learning to Optimize: A Primer and A Benchmark. Journal of Machine Learning Research, 23(189):1–59, 2022. [19] X. Chen, T. Chen, Y. Cheng, W. Chen, A. Awadallah, and Z. Wang. Scalable Learning to Optimize: A Learned Optimizer Can Train Big Models. In Computer Vision – ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXIII, pages 389–405, Berlin, Heidelberg, Oct. 2022. Springer-Verlag. [20] X. Chen, J. Liu, Z. Wang, and W. Yin. Theoretical Linear Convergence of Unfolded ISTA and its Practical Weights and Thresholds. In Advances in Neural Information Processing Systems, volume 31, pages 9079–9089, Canada, 2018. Curran Associates, Inc. 15
[21] S. Das Gupta, B. P. G. Van Parys, and E. K. Ryu. Branch-and-bound performance estimation programming: A unified methodology for constructing optimal optimization methods. Mathematical Programming, 204(1-2):567–639, Mar. 2024. [22] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei. Imagenet: A large-scale hierarchical image database. In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. IEEE, 2009. [23] S. Diamond and S. Boyd. Total variation inpainting. CVXPY Examples, 2016. Accessed: 2025-05-01. [24] Y. Drori and M. Teboulle. Performance of first-order methods for smooth convex minimization: A novel approach. Mathematical Programming, 145(1):451–482, June 2014. [25] G. K. Dziugaite and D. M. Roy. Computing Nonvacuous Generalization Bounds for Deep (Stochastic) Neural Networks with Many More Parameters than Training Data. In Proceedings of the 33rd Conference on Uncertainty in Artificial Intelligence. AUAI Press, Aug. 2017. [26] N. Fournier and A. Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162(3-4):707–738, Aug. 2015. [27] M. Frank and P. Wolfe. An algorithm for quadratic programming. Naval Research Logistics Quarterly, 3(1-2):95–110, 1956. [28] M. Garstka, M. Cannon, and P. Goulart. COSMO: A conic operator splitting method for large convex problems. In 2019 18th European Control Conference (ECC), pages 1951–1956, Naples, Italy, June 2019. IEEE. [29] J.-y. Gotoh, M. J. Kim, and A. E. B. Lim. Calibration of Distributionally Robust Empirical Optimization Models. Operations Research, 69(5):1630–1650, Sept. 2021. [30] P. J. Goulart and Y. Chen. Clarabel: An interior-point solver for conic programs with quadratic objectives, May 2024. [31] K. Gregor and Y. LeCun. Learning fast approximations of sparse coding. In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, pages 399–406, Madison, WI, USA, June 2010. Omnipress. [32] Q. Healey, P. Nobel, and S. Boyd. Differentiating Through a Quadratic Cone Program, Aug. 2025. [33] H. Heaton, X. Chen, Z. Wang, and W. Yin. Safeguarded Learned Convex Optimization. Proceedings of the AAAI Conference on Artificial Intelligence, 37(6):7848–7855, June 2023.
16
[34] Z. Hu and L. J. Hong. Kullback-Leibler Divergence Constrained Distributionally Robust Optimization, Nov. 2012. [35] J. Huang, P. Goulart, and K. Margellos. Data-Driven Performance Guarantees for Parametric Optimization Problems, June 2025. [36] U. Jang, S. D. Gupta, and E. K. Ryu. Computer-Assisted Design of Accelerated Composite Optimization Methods: OptISTA, Sept. 2024. [37] Y. Kamri, J. M. Hendrickx, and F. Glineur. Numerical Design of Optimized First-Order Algorithms, July 2025. [38] D. Kim. Accelerated proximal point method for maximally monotone operators. Mathematical Programming, 190(1):57–87, Nov. 2021. [39] D. Kim and J. A. Fessler. Optimized first-order methods for smooth convex minimization. Mathematical Programming, 159(1-2):81–107, Sept. 2016. [40] D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. In 3rd International Conference on Learning Representations, 2015. [41] D. Kuhn, P. M. Esfahani, V. A. Nguyen, and S. Shafieezadeh-Abadeh. Wasserstein Distributionally Robust Optimization: Theory and Applications in Machine Learning. In S. Netessine, D. Shier, and H. J. Greenberg, editors, Operations Research & Management Science in the Age of Analytics, pages 130–166. INFORMS, Washington, Oct. 2019. [42] D. Kuhn, S. Shafiee, and W. Wiesemann. Distributionally Robust Optimization. Acta Numerica, 34:579–804, 2025. [43] Y. Le and X. Yang. Tiny ImageNet visual recognition challenge. CS 231N course report, Stanford University, 2015. [44] L. Lessard, B. Recht, and A. Packard. Analysis and Design of Optimization Algorithms via Integral Quadratic Constraints. SIAM Journal on Optimization, 26(1):57–95, Jan. 2016. [45] K. Li and J. Malik. Learning to Optimize. In International Conference on Learning Representations, Feb. 2017. [46] Z. Lin, Z. Xiong, D. Ge, and Y. Ye. A Practical GPU-Enhanced Matrix-Free PrimalDual Method for Large-Scale Conic Programs, Apr. 2026. [47] J. Liu, X. Chen, Z. Wang, and W. Yin. ALISTA: Analytic Weights Are As Good As Learned Weights in LISTA. In International Conference on Learning Representations, May 2019.
17
[48] J. Liu, X. Chen, Z. Wang, W. Yin, and H. Cai. Towards Constituting Mathematical Structures for Learning to Optimize. In Proceedings of the 40th International Conference on Machine Learning, pages 21426–21449. PMLR, July 2023. [49] I. Loshchilov and F. Hutter. SGDR: Stochastic Gradient Descent with Warm Restarts. In International Conference on Learning Representations, 2017. [50] I. Loshchilov and F. Hutter. Decoupled Weight Decay Regularization. In International Conference on Learning Representations, May 2019. [51] H. Lu, Z. Peng, and J. Yang. cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear Programming, Sept. 2025. [52] H. Lu and J. Yang. cuPDLP.jl: A GPU Implementation of Restarted Primal-Dual Hybrid Gradient for Linear Programming in Julia, June 2024. [53] V. A. Marčenko and L. A. Pastur. Distribution of Eigenvalues for Some Sets of Random Matrices. Mathematics of the USSR-Sbornik, 1(4):457–483, Apr. 1967. [54] A. Martin and G. Belgioioso. Learning to accelerate Krasnosel’skii-Mann fixed-point iterations with guarantees, Jan. 2026. [55] A. Martin and L. Furieri. Learning to Optimize With Convergence Guarantees Using Nonlinear System Theory. IEEE Control Systems Letters, 8:1355–1360, 2024. [56] A. Martin, I. R. Manchester, and L. Furieri. Learning to optimize with guarantees: A complete characterization of linearly convergent algorithms, Aug. 2025. [57] P. Mohajerin Esfahani and D. Kuhn. Data-driven distributionally robust optimization using the Wasserstein metric: Performance guarantees and tractable reformulations. Mathematical Programming, 171(1-2):115–166, Sept. 2018. [58] B. O’Donoghue. Operator splitting for a homogeneous embedding of the linear complementarity problem. SIAM Journal on Optimization, 31(3):1999–2023, 2021. [59] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd. Conic optimization via operator splitting and homogeneous self-dual embedding. Journal of Optimization Theory and Applications, 169(3):1042–1068, 2016. [60] N. Parikh and S. Boyd. Proximal Algorithms. Foundations and Trends® in Optimization, 1(3):127–239, 2014. [61] J. Park, V. Ranjan, and B. Stellato. Data-driven Analysis of First-Order Methods via Distributionally Robust Optimization, Dec. 2025. [62] F. Pedregosa and D. Scieur. Average-case Acceleration through spectral density estimation. In Proceedings of the 37th International Conference on Machine Learning, pages 7553–7562, Austria (Virtual), Nov. 2020. PMLR. 18
[63] I. Prémont-Schwarz, J. Vı́tků, and J. Feyereisl. A Simple Guard for Learned Optimizers. In Proceedings of the 39th International Conference on Machine Learning, pages 17910– 17925. PMLR, July 2022. [64] V. Ranjan, J. Park, S. Gualandi, A. Lodi, and B. Stellato. Exact Verification of FirstOrder Methods via Mixed-Integer Linear Programming, May 2025. [65] V. Ranjan and B. Stellato. Verification of first-order methods for parametric quadratic optimization. Mathematical Programming, July 2025. [66] R. T. Rockafellar and S. Uryasev. Conditional Value-at-Risk: Optimization Approach. In S. Uryasev and P. M. Pardalos, editors, Stochastic Optimization: Algorithms and Applications, pages 411–435. Springer US, Boston, MA, 2001. [67] E. K. Ryu, A. B. Taylor, C. Bergeling, and P. Giselsson. Operator Splitting Performance Estimation: Tight Contraction Factors and Optimal Parameter Selection. SIAM Journal on Optimization, 30(3):2251–2271, Jan. 2020. [68] E. K. Ryu and W. Yin. Large-Scale Convex Optimization: Algorithms & Analyses via Monotone Operators. Cambridge University Press, Cambridge, UK, 1 edition, 2022. [69] F. Samaria and A. Harter. Parameterisation of a stochastic model for human face identification. In Proceedings of 1994 IEEE Workshop on Applications of Computer Vision, pages 138–142, 1994. [70] R. Sambharya, J. Bok, N. Matni, and G. Pappas. Learning Acceleration Algorithms for Fast Parametric Convex Optimization with Certified Robustness, Oct. 2025. [71] R. Sambharya and B. Stellato. Learning Algorithm Hyperparameters for Fast Parametric Convex Optimization, Nov. 2024. [72] R. Sambharya and B. Stellato. Data-Driven Performance Guarantees for Classical and Learned Optimizers. Journal of Machine Learning Research, 26(171):1–49, 2025. [73] S. Smale. On the average number of steps of the simplex method of linear programming. Mathematical Programming, 27(3):241–262, Oct. 1983. [74] Q. Song, W. Lin, J. Wang, and H. Xu. Towards Robust Learning to Optimize with Theoretical Guarantees. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 27498–27506, 2024. [75] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd. OSQP: An operator splitting solver for quadratic programs. Mathematical Programming Computation, 12(4):637–672, 2020. [76] M. Sucker, J. Fadili, and P. Ochs. Learning-to-Optimize with PAC-Bayesian Guarantees: Theoretical Considerations and Practical Implementation. Journal of Machine Learning Research, 26(211):1–53, 2025. 19
[77] A. B. Taylor, J. M. Hendrickx, and F. Glineur. Smooth strongly convex interpolation and exact worst-case performance of first-order methods. Mathematical Programming, 161(1-2):307–345, Jan. 2017. [78] R. Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2018. [79] W. Wiesemann, D. Kuhn, and M. Sim. Distributionally Robust Convex Optimization. Operations Research, 62(6):1358–1376, Dec. 2014. [80] Z. Xiong. High-Probability Polynomial-Time Complexity of Restarted PDHG for Linear Programming, Jan. 2025. [81] J. Yang, X. Chen, T. Chen, Z. Wang, and Y. Liang. M-L2O: Towards Generalizable Learning-to-Optimize by Test-Time Fast Self-Adaptation. In The Eleventh International Conference on Learning Representations, May 2023.
20
A
Additional details and proofs for Section 3
A.1
Tractable convex conic formulation of (PEP)
The performance estimation problem (PEP) admits the following tractable formulation: sup ℓ(Aθ (z)) = maximize ℓ(G, F ) = tr(ATobj G) + bTobj F z∈Z subject to S(G, F ) ∈ RM + tr(AT0 G) + bT0 F + c0 ≤ 0 G ∈ SK+2 , F ∈ RK+1 . + The Lagrangian dual of the problem above writes as minimize −c0 τ subject to −S ∗ (y) + τ (A0 , b0 ) − (Aobj , bobj ) ∈ SK+2 × {0} + τ ≥ 0, y ∈ RM +, ∗ M where → SK+2 × RK+1 is an adjoint mapping of S defined as S ∗ (y) = PM S : R − m=1 ym (Am , bm ). For the most of the first-order methods, strong duality holds true from [77, Theorem 6].
A.2
Tractable form of (DRO-PEP) with expected loss
The problem (DRO-PEP) with expected loss admits a tractable convex conic formulation [61]: P minimize (1/N ) N i=1 si bi , Fbi ) + λε ≤ si subject to −c0 τi − (Xi , Yi ), (G −S ∗ (yi ) − (Xi , Yi ) + τi (A0 , b0 ) − (Aobj , bobj ) ∈ SK+2 × {0} + ∥(Xi , Yi )∥∗ ≤ λ, i = 1, . . . , N,
(DRO-PEP-D)
with optimization variables si ∈ R, τi ∈ R+ , (Xi , Yi ) ∈ SK+2 × RK+1 , yi ∈ RM + for i = 1, . . . , N , and λ ∈ R+ , where S ∗ is an adjoint mapping of S and ∥ · ∥∗ is a dual norm of ∥ · ∥ over SK+2 × RK+1 .
A.3
Proof of Proposition 1
We start the proof with the strong duality result between (DRO-PEP) and its dual (DROPEP-D). Our proof reformulates that of [57, Corollary 5.1]. Lemma 1. Problem (DRO-PEP-D) is a dual of (DRO-PEP) and strong duality holds whenever ε > 0. 21
Proof. The statement of the lemma is exactly the same as that of [57, Corollary 5.1], except for showing the strong duality result for the following problem: σΞθ (ν) = sup ν, ξ s.t. ξ ∈ Ξθ . According to the strong duality result of [77, Theorem 6] where θ ∈ Θ satisfies the assumption that such a step size only represents algorithms that uses information of g k for xk+1 -update, the problem above has zero duality gap and the evaluation is exact with its (feasible) dual problem formulation. ■ Next, we show that the bi-dual of (DRO-PEP), which is a dual form for (DRO-PEP-D), is exactly the form of (DRO-PEP-P) and strong duality holds. Lemma 2. Problem (DRO-PEP-P) is a dual of (DRO-PEP-D) and strong duality holds whenever ε > 0. Proof. The fact that (DRO-PEP-P) is a dual of (DRO-PEP-D) comes from [61, Theorem 2]. The strong duality result comes from the fact that (DRO-PEP-D) admits a slater point for any ε > 0, which is exactly the result of [61, Lemma 4]. ■ Proof of Proposition 1. Combining Lemma 1 and Lemma 2 gives the desired result.
A.4
■
Extended dual formulation of (DRO-PEP) with linear operator interpolations
Consider the case where the support set Ξ is defined with additional positive semidefinite constraint as follows. T T tr(A G) + b F + c ≤ 0 0 0 0 K+2 K+1 T T Ξ = (G, F ) ∈ S+ × R . tr(Am G) + bm F ≤ 0, m = 1, . . . , M, K+1 T T H ∈ S+ where (H)kl = tr(Ckl G) + dkl F The class of linear operators with bounded spectral norm and convex quadratic functions are such cases, as the following example illustrates. Example 1 (Interpolation condition of linear operators). We state [12, Theorem 3.1], which provides with interpolation condition for linear operators, along with convex quadratic function as its special case. Let X ∈ Rn×N1 , Y ∈ Rm×N1 , U ∈ Rm×N2 , V ∈ Rn×N2 and L ≥ 0. Then there exists M ∈ Rn×m with σmax (M ) ≤ L such that U = M T Y,
V = M X, if and only if X T V = Y T U,
Y T Y ⪯ L2 X T X, 22
V T V ⪯ L2 U T U.
Note that the first condition can be written as tr(ATm G)+bTm F ≤ 0 and − tr(ATm G)−bTm F ≤ 0 by writing each entries of the matrix X T P (or P T X) as the left-hand side of these inequalities. For the second constraint, let H = L2 X T X − Y T Y and H2 = L2 U T U − V T V . Then H ⪰ 0 and each entries of H are linear combinations of bilinear terms of (X, Y ). Therefore, it is of the form T (H)kl = tr(Ckl G) + dTkl F, 1 ≤ k, l ≤ K + 1,
where Ckl ∈ SK+2 and dkl ∈ RK+1 . We can do similarly for the last constraint. The conditions above can be specialized to symmetric positive semidefinite matrix Q such that that µI ⪯ Q ⪯ LI, in a sense that g k = Qxk and f k = (1/2)(xk )T Qxk for k = 1, . . . , K if any only if, for X = (x0 , . . . , xK ) ∈ Rd×(K+1) , P = (g 0 , . . . , g K ) ∈ Rd×(K+1) , and F = (f 0 , . . . , f K ) ∈ RK+1 , (X, P, F ) is an element of following set: XT P = P T X d×(K+1) d×(K+1) K+1 T Qµ,L = (X, P, F ) ∈ R ×R ×R . (P − µX) (LX − P ) ⪰ 0 T fi = (1/2)xi gi , i = 0, . . . , K (3)
Lemma 3. Consider the class of convex quadratic functions of (3) as our choice of F. The worst-case expected loss of the form maximize
E ℓ(Aθ (z))
z∼Q
b N ), subject to Q ∈ Uε (P
supp Q = Ξ,
is the optimal value of the problem P minimize (1/N ) N i=1 si bi ) − YiT Fbi + λε ≤ si subject to −c0 τi − tr(XiT G PK+1 e P K+2 τi A0 + M k,l=1 (Hi )kl Ckl − Aobj − Xi ∈ S+ m=1 (yi )m Am − PK+1 PM e i )kl dkl − bobj − Yi = 0 τi b0 + m=1 (yi )m bm − k,l=1 (H ∥(Xi , Yi )∥ ≤ λ, i = 1, . . . , N, K+2 e i ∈ SK+1 with variables si ∈ R, τi ∈ R+ , yi ∈ RM , Yi ∈ RK+1 , H , and λ ∈ R+ . + + , Xi ∈ S
Proof. The worst-case expected loss given can be written as sup E ℓ(Aθ (z)) = sup E bN) Q∈Uε (P
z∼Q
bθ ) Qθ ∈Uε (P N
(G,F )∼Qθ
ℓ(G, F ) .
Circling back to the proof of [61, Theorem 2], this is equivalent to P minimize (1/N ) N i=1 si ∗ bi ) − YiT Fbi + λε ≤ si subject to (−ℓ) (Xi , Yi ) − (Ui , Vi ) + σΞ (Ui , Vi ) − tr(XiT G ∥(Xi , Yi )∥ ≤ λ, i = 1, . . . , N, 23
with variables si ∈ R, λ ∈ R+ , Xi , Ui ∈ SK+2 , and Yi , Vi ∈ RK+1 for i = 1, . . . , N . First of all, (−ℓ)∗ (X, Y ) = sup tr(X T G) + Y T F + ℓ(G, F ) (G,F )∈SK+2 ×RK+1
=
sup
tr (X + Aobj )T G + (Y + bobj )T F
(G,F )∈SK+2 ×RK+1
( 0 if (X, Y ) + (Aobj , bobj ) = 0 = ∞ otherwise. Furthermore, the support function σΞ of Ξ is σΞ (U, V ) = sup tr(U T G) + V T F s.t. (G, F ) ∈ Ξ
= sup tr(U T G) + V T F s.t. tr(AT0 G) + bT0 F + c0 ≤ 0 tr(ATm G) + bTm F ≤ 0, m = 1, . . . , M T H ⪰ 0, (H)kl = tr(Ckl G) + dTkl F (G, F ) ∈ SK+2 × RK+1 +
≤ inf −c0 τ PK+1 P K+2 e s.t. τ A0 + M k,l=1 Ckl Hkl − U ∈ S+ m=1 ym Am − PK+1 PM e kl − V = 0, τ b0 + m=1 ym bm − k,l=1 dkl H K+1 e with variables τ ∈ R+ , y ∈ RM . + , and H ∈ S+ According to the example in [77, Theorem 6], there exists a full-rank G and F such that the Slater’s condition holds for PEP with class of L-smooth µ-strongly convex functions. Using the G ≻ 0 and F build for (µ + ε)-strongly convex (L − ε)-smooth functions for some ε > 0 such that µ + ε ≤ L − ε, this is actually G ≻ 0 and the positive semidefinite constraint of (3) strictly holds. According to the Slater’s condition, the inequality above is actually an equality. Combining the results above, the worst-case expected loss is the optimal value of the problem P minimize (1/N ) N i=1 si bi ) − YiT Fbi + λε ≤ si subject to −c0 τi − tr(XiT G P PK+1 e K+2 τi A0 + M m=1 (yi )m Am − k,l=1 (Hi )kl Ckl − Aobj − Xi ∈ S+ PM PK+1 e i )kl dkl − bobj − Yi = 0 τi b0 + m=1 (yi )m bm − k,l=1 (H ∥(Xi , Yi )∥ ≤ λ, i = 1, . . . , N, K+2 e i ∈ SK+1 with variables si ∈ R, τi ∈ R+ , yi ∈ RM , Yi ∈ RK+1 , H , and λ ∈ R+ + + , Xi ∈ S and the strong duality holds. ■
24
B
Proofs for Section 5
B.1
Proof of Theorem 1
First of all, we show the Lipschitz property of algorithm mapping θ 7→ Aθ . Lemma 4. Let K and L > 0 be fixed. Suppose Assumption 2 holds, the algorithm Aθ involves either gradient step of L-smooth function f : x
k+1
k
=x +
k X i=0
θk,i ∇f (xi ),
or proximal step of convex, proper, and lower semi-continuous g: xk+1 = proxθg (xk ), and ∇f (x⋆ ) and x⋆ − proxηg x⋆ are uniformly bounded over choices of f and g. Then the mapping θ 7→ Aθ is Lipschitz, i.e., there exists L̄ > 0 such that ∥Aθ1 (z) − Aθ2 (z)∥ ≤ L̄ ∥θ1 − θ2 ∥ , for all θ1 , θ2 ∈ Θ and z ∈ Z.
T 0 Proof. Suppose that the grammian representation G ∈ SK + is defined as G = P P where P ∈ Rd×K0 is a horizontal stack of vectors of the forms: xi −x⋆ , ∇f (xi ), or proxθg (xi ). We denote by G = Aθ (z) = G(θ), P = P (θ), and xi = xi (θ) to emphasize the θ-dependency of each terms. For the sake of generality, let x⋆ = 0 for all θ ∈ Θ. Since
Aθ1 (z) − Aθ2 (z) = P1T P1 − P2T P2 =
1 (P1 − P2 )T (P1 + P2 ) + (P1 + P2 )T (P1 − P2 ) , 2
we have ∥Aθ1 (z) − Aθ2 (z)∥F ≤ ∥P (θ1 ) − P (θ2 )∥F ∥P (θ1 ) + P (θ2 )∥2 n o ≤ 2 ∥P (θ1 ) − P (θ2 )∥F max max xi (θ) − x⋆ , ∇f (xi (θ)), proxθg (xi (θ)) . θ∈{θ1 ,θ2 }
i
It remains to show that θ 7→ P (θ) is Lipschitz continuous and the vectors xi (θ) − x⋆ , ∇f (xi (θ)), and proxθg (xi (θ))}K i=0 are uniformly bounded for i = 0, . . . , K. We first prove the latter claim by the following recursive statement: If xk (θ) − x⋆ is bounded, then xk+1 (θ) − x⋆ is bounded for these two cases: xk+1 (θ) = xk (θ) − P k k k+1 (θ) = proxθk g (xk (θ)). Note that we impose the initial condii=0 θk,i ∇f (x (θ)) or x tion that x0 (θ) − x⋆ is uniformly bounded with θ ∈ Θ. Also, from L-Lipschitz of ∇f , ∇f (xk (θ)) − ∇f (x⋆ ) ≤ L xk (θ) − x⋆ so ∇f (xk (θ)) is bounded as well. For the first case, we have xk+1 (θ) − x⋆ 25
k
⋆
x (θ) − x
=
k
⋆
≤ x (θ) − x
−
k X
+
i=0
⋆
θk,i ∇f (x (θ)) − ∇f (x ) −
k X i=0
i
k X i=0
θk,i ∇f (x⋆ )
k X θk,i ∇f (x (θ)) − ∇f (x ) + θk,i ∥∇f (x⋆ )∥ i
⋆
i=0
≤ xk (θ) − x⋆
v v u k u k k uX uX X 2 t 2 t i ⋆ θ + θ ∥∇f (x (θ)) − ∇f (x )∥ +
≤ xk (θ) − x⋆
v u k k uX X 2 t i ⋆ + diam(Θ) θk,i ∥∇f (x⋆ )∥ ∥∇f (x (θ)) − ∇f (x )∥ +
≤ xk (θ) − x⋆
v u k k uX X θk,i ∥∇f (x⋆ )∥ . + L diam(Θ)t ∥xi (θ) − x⋆ ∥2 +
k,i
k,i
i=0
i=0
i=0
∥∇f (x⋆ )∥
i=0
i=0
i=0
i=0
For the second case, xk+1 (θ) − x⋆ =
proxθk g (xk (θ)) − proxθk g (x⋆ ) + proxθk g (x⋆ ) − x⋆
≤ proxθk g (xk (θ)) − proxθk g (x⋆ ) + proxθk g (x⋆ ) − x⋆ ≤ xk (θ) − x⋆ + proxθk g (x⋆ ) − x⋆ ,
where the inequality comes from proxθk g being a firmly-nonexpansive operator. Therefore, all iterates are bounded uniformly over θ ∈ Θ, given a fixed K. Since G = P T P with P S stacking these bounded iterates and gradients, the support set ΞΘ = θ∈Θ Ξθ is bounded. Now it remains to show the Lipschitz property of θ 7→ P (θ). We similarly use induction to prove that, if xk (θ1 ) − xk (θ2 ) ≤ Lk ∥θ1 − θ2 ∥ for all θ1 , θ2 ∈ Θ and any choice of f and g, then there exists Lk+1 > 0 such that xk+1 (θ1 ) − xk+1 (θ2 ) ≤ Lk+1 ∥θ1 − θ2 ∥ for all θ1 , θ2 ∈ Θ and any f and g as well. For the first case (f -related update), we have xk+1 (θ1 ) − xk+1 (θ2 ) = xk (θ1 ) − xk (θ2 ) − = xk (θ1 ) − xk (θ2 ) − = xk (θ1 ) − xk (θ2 ) − +
k X i=0
k X i=0 k X i=0 k X i=0
(θ1 )k,i ∇f (xi (θ1 )) − (θ2 )k,i ∇f (xi (θ2 ))
k X (θ1 )k,i − (θ2 )k,i ∇f (xi (θ1 )) − (θ2 )k,i ∇f (xi (θ1 )) − ∇f (xi (θ2 )) i=0
(θ1 )k,i − (θ2 )k,i ∇f (xi (θ1 )) − ∇f (x⋆ )
(θ1 )k,i − (θ2 )k,i
⋆
∇f (x ) −
k X i=0
(θ2 )k,i ∇f (xi (θ1 )) − ∇f (xi (θ2 )) 26
Therefore, xk+1 (θ1 ) − xk+1 (θ2 ) ≤ xk (θ1 ) − xk (θ2 ) + + +
k X
i=0 k X i=0 k X i=0
(θ1 )k,i − (θ2 )k,i ∇f (xi (θ1 )) − ∇f (x⋆ ) (θ1 )k,i − (θ2 )k,i ∇f (x⋆ )
(θ2 )k,i ∇f (xi (θ1 )) − ∇f (xi (θ2 ))
k
≤ x (θ1 ) − xk (θ2 ) v v u k u k uX uX 2 t (θ ) − (θ ) ∥∇f (xi (θ )) − ∇f (x⋆ )∥2 +t 1 k,i
2 k,i
1
i=0
i=0
v u k uX +t (θ ) i=0
1 k,i − (θ2 )k,i
2
∥∇f (x⋆ )∥
v v u k u k uX uX ∥∇f (xi (θ )) − ∇f (xi (θ ))∥2 + t (θ )2 t 2 k,i
i=0
1
2
i=0
v u k uX k k ≤ x (θ1 ) − x (θ2 ) + L ∥θ1 − θ2 ∥ t ∥xi (θ1 ) − x⋆ ∥2 i=0
v u k uX ⋆ + ∥θ1 − θ2 ∥∥∇f (x )∥ + L diam(Θ)t ∥xi (θ1 ) − xi (θ2 )∥2 . i=0
For the second case (g-related update), we have xk+1 (θ1 ) − xk+1 (θ2 )
= prox(θ1 )k g (xk (θ1 )) − prox(θ2 )k g (xk (θ2 ))
= prox(θ1 )k g (xk (θ1 )) − prox(θ1 )k g (xk (θ2 )) + prox(θ1 )k g (xk (θ2 )) − prox(θ2 )k g (xk (θ2 ))
≤ prox(θ1 )k g (xk (θ1 )) − prox(θ1 )k g (xk (θ2 )) + prox(θ1 )k g (xk (θ2 )) − prox(θ2 )k g (xk (θ2 )) ≤ xk (θ1 ) − xk (θ2 ) + prox(θ1 )k g (xk (θ2 )) − prox(θ2 )k g (xk (θ2 )) .
If g is an indicator of closed convex set, therefore proxθg is always a projection mapping, then xk+1 (θ1 ) − xk+1 (θ2 ) ≤ xk (θ1 ) − xk (θ2 ) . 27
For other cases, we have xk (θ2 ) − prox(θ2 )k g (xk (θ2 )) prox(θ1 )k g (x (θ2 )) − prox(θ2 )k g (x (θ2 )) ≤ |(θ1 )k − (θ2 )k | . (θ2 )k k
k
Note that for any α > 0, x − proxαg (x) ∈ ∂g proxαg (x) . α We have already shown that the iterates are uniformly bounded over θ ∈ Θ and the choice of g. Therefore, the iterates lie on compact set of Rd . As we can restrict the domain of g to a compact set, ∂g(proxαg (x)) for any x within such compact set is always locally bounded [8, Proposition 16.17]. We can choose finite covering of Θ in terms of such neighborhoods, which translates to a finite covering of compactly-restricted domain of g. Within such compactlyrestricted domain, ∂g is uniformly bounded over θ ∈ Θ, which results in sup θ∈Θ
xk (θ) − proxθk g (xk (θ)) ≤ M < ∞. θk
Therefore, xk+1 (θ1 ) − xk+1 (θ2 ) ≤ xk (θ1 ) − xk (θ2 ) + ∥θ1 − θ2 ∥ R,
for some universal constant R > 0 over θ1 , θ2 ∈ Θ and g. We may conclude that there exists Lk+1 > 0 such that, for all θ1 , θ2 ∈ Θ and functions f , g, xk+1 (θ1 ) − xk+1 (θ2 ) ≤ Lk+1 ∥θ1 − θ2 ∥ , given xk (θ1 ) − xk (θ2 ) P for L̄ = K k=0 Lk .
≤ Lk ∥θ1 − θ2 ∥. This results in ∥P (θ1 ) − P (θ2 )∥F ≤ L̄ ∥θ1 − θ2 ∥ ■
Remark 2 (Assumption in Lemma 4). We have seemingly nontrivial assumption on Lemma 4 that, given η > 0, ∇f (x⋆ ) and x⋆ − proxηg x⋆ being uniformly bounded over choices of f and g. This can be implied by the condition that there exist x⋆f ∈ argminx f (x) and x⋆g ∈ argminx f (x) such that x⋆f and x⋆g are not too far away from x⋆ , as ∥∇f (x⋆ )∥ = ∇f (x⋆ ) − ∇f (x⋆f ) ≤ L x⋆ − x⋆f , and x⋆ − proxηg (x⋆ ) = (x⋆ − proxηg (x⋆ )) − (x⋆g − proxηg (x⋆g )) ≤ x⋆ − x⋆g . We follow the measure concentration result of [26] to derive the radius ε guaranteeing that DRO loss upper-bounds true risk with probability at least 1−β, as in [61]. Choose ε = ε(β, θ) as the follows: ((K+2)(K+3)/2+(K+1))−1 log(c1 β −1 ) log(c1 β −1 ) N ≥ c2 ε(β, θ) = c2 N−1 1/a (4) −1 ) log(c β ) log(c 1 1β N< , c2 N c2 28
where c1 > 0 and c2 > 0 are some appropriately chosen constants depending only on θ, due to its dependency on true distribution Pθ , which is the pushforward measure of P by mapping Aθ . From this, we get the following measure concentration inequality for each θ ∈ Θ: Lemma 5. Let β ∈ (0, 1). For each θ ∈ Θ, there exists ε(β, θ) > 0 as in (4) such that bθ ≥ 1 − β. PN Pθ ∈ Uε(β,θ) P N Proof. Proof follows directly from [26, Theorem 2].
■
We now prove Theorem 1 using the covering number argument on compact Θ to extract θ-independent radius ε of Wasserstein ambiguity set. Proof of Theorem 1. For any δ > 0, there exists a finite covering Cδ of compact set Θ in a sense that, |Cδ | < ∞ and for any θ ∈ Θ, there exists θδ ∈ Cδ such that ∥θ − θδ ∥ ≤ δ. According to Lemma 5, there exists ε(β/|Cδ |) = supθj ∈Cδ ε(β/|Cδ |, θj ) > 0 such that P
N
b θj P ∈ Uε(β/|Cδ |) P N θj
≥ 1 − β/|Cδ |,
∀θj ∈ Cδ .
From the union bound, we have b θj , ∀θj ∈ Cδ ≥ 1 − β. PN Pθj ∈ Uε(β/|Cδ |) P N Now, define ε(β, δ) = ε(β/|Cδ |) + 2Lδ. For any θ ∈ Θ, there exists θj ∈ Cδ such that ∥θ − θj ∥ ≤ δ. From Lemma 4, we have bθ , P b θj ≤ L ∥θ − θj ∥ ≤ Lδ, W1 P N N and W1 Pθ , Pθj ≤ L ∥θ − θj ∥ ≤ Lδ. From triangle inequality of Wasserstein distance, b θ , Pθ W1 P N θj θj θ θj b b b ≤ W1 PN , PN + W1 PN , P + W 1 P θj , P θ ≤ Lδ + ε(β/|Cδ |) + Lδ = ε(β, δ). Therefore, we get b θ ), ∀θ ∈ Θ ≥ 1 − β. PN Pθ ∈ Uε(β,δ) (P N 29
Finally, for any δ > 0, we get θj θj θ N N θ b b P P ∈ Uε(β) PN , ∀θ ∈ Θ ≥ P P ∈ Uε(β/|Cδ |) PN , ∀θj ∈ Cδ ≥ 1 − β, for ε(β, δ) = ε(β/|Cδ |) + 2Lδ. From [78, Corollary 4.2.13], the covering number |Cδ | of Θ is upper bounded as |Cδ | ≤
2 diam(Θ) +1 δ
dim(Θ) .
Then ε(β, δ) ≥ ε β(2 diam(Θ)/δ + 1)−dim(Θ) + 2Lδ, so the choice of δ with tightest ε(β, δ) will give the best bound on ε(β). From (4), we have ε(β, δ) ≥
log(1/β) + dim(Θ) log(2 diam(Θ)/δ + 1) N
1/q + 2Lδ,
where q is the dimension of (G, F ). Naively choosing δ = 1/N yields ε(β, 1/N ) ≲
log(1/β) + dim(Θ) log diam(Θ) + log N N
!1/q +
2L . N ■
B.2
Proof of Proposition 2
Throughout this subsection, we omit θ-dependency and write Pθ and Qθ as P and Q, for simplicity. First, we state the proposition from [42], which will be useful in our proofs. Proposition 3 (Proposition 8.5 of [42]). Consider the 1-Wasserstein ambiguity set n o b N ) = Q | W1 (Q, P bN) ≤ ε , Uε (P where Q is a probability measure supported on a closed set Ξ. Then, if ℓ is a Lipschitz continuous function with EY ∼Pb N (|ℓ(Y )|) < ∞, sup bN) Q∈Uε (P
E ℓ(Y ) ≤
Y ∼Q
E
ℓ(Y ) + ε Lip(ℓ).
bN Y ∼P
Note that from [61], DRO loss interpolates (OPT-PEP) and (L2O) in the following sense: the problem (DR-L2O) either converges to (OPT-PEP) when ε → ∞ or is equivalent to (L2O) when ε = 0, given a fixed θ ∈ Θ. This property is directly stated in the following lemma. 30
Lemma 6. For any algorithm parameter θ ∈ Θ and probability distribution Q supported on Ξ, R(θ, Q) ≤ Rε (θ, Q) ≤ Rε+ (θ, Q) ≤ sup ℓ(G, F ), (G,F )∈Ξ
for any ε+ ≥ ε ≥ 0. Proof. From {Q} ⊆ Uε (Q) ⊆ Uε+ (Q), we get R(θ, Q) ≤ Rε (θ, Q) ≤ Rε+ (θ, Q). If Ξ is compact, a maximizer (G⋆ , F ⋆ ) ∈ argmax(G,F )∈Ξ ℓ(G, F ) exists. Therefore, for every ε > 0, ℓ(G, F ) ≤ sup sup ℓ(G, F ) = sup ℓ(G, F ). Rε (θ, Q) = sup E e (G,F )∼Q e Q∈U ε (Q)
e Q∈U ε (Q) (G,F )∈Ξ
(G,F )∈Ξ
■ We now prove Proposition 2. Proof of Proposition 2. First, from Lemma 6, b N ) ≤ Rε (θ, P b N ) ≤ Rε+ (θ, P b N ) ≤ sup ℓ Aθ (z) , R(θ, P z∈Z
for any θ ∈ Θ and ε+ ≥ ε > 0. Since Θ is a compact set, each terms in inequality above has respective minimizers in Θ. First of all, b N ) ≤ R(θ, P b N ) ≤ Rε (θ, P b N ) ≤ Rε+ (θ, P b N ), inf R(θ, P
θ∈Θ
∀θ ∈ Θ.
b N ), we have For θε⋆ ∈ argminθ∈Θ Rε (θ, P b N ) ≤ R(θ⋆ , P b N ) ≤ Rε (θ⋆ , P b N ) = inf Rε (θ, P b N ) ≤ Rε (θ, P b N ), inf R(θ, P ε ε
θ∈Θ
θ∈Θ
∀θ ∈ Θ.
b N ), Similarly, with respect to θε+ ∈ argminθ∈Θ Rε+ (θ, P b N ) ≤ Rε (θ⋆ , P b N ) ≤ Rε+ (θ⋆+ , P b N ) ≤ Rε+ (θ, P b N ), inf R(θ, P ε ε
θ∈Θ
∀θ ∈ Θ.
Lastly, b N ) ≤ Rε (θ⋆ , P b N ) ≤ Rε+ (θ⋆+ , P b N ) ≤ inf inf R(θ, P ε ε
θ∈Θ
sup
θ∈Θ (G,F )∈Ξθ
ℓ(G, F ),
b N ) is monotonically nondecreasing and upper-bounded we observe that ε 7→ Rε (θε⋆ , P by inf θ∈Θ supz∈Z ℓ Aθ (z) . b N ) ≤ Rε (θ, P b N ) implies inf θ∈Θ R(θ, P b N ) ≤ inf θ∈Θ Rε (θ, P b N ), the Furthermore, as R(θ, P bN) ≤ following holds true for all ε > 0. The left inequality follows from Lemma 6 (since R(θ, P 31
b N ) for all θ, taking inf θ∈Θ on both sides gives inf θ R(θ, P b N ) ≤ inf θ Rε (θ, P bN) = Rε (θ, P ⋆ b Rε (θε , PN )). The right inequality follows from Proposition 3 (applying Proposition 8.5 of [42] to each θ and taking inf θ∈Θ on the right-hand side): b N ) ≤ Rε (θ⋆ , P b N ) ≤ inf R(θ, P b N ) + ε Lip(ℓ). inf R(θ, P ε
θ∈Θ
θ∈Θ
Applying ε → 0+ , we get the desired result. b N ) is (jointly) continuous in (ε, θ). Also, Note that mapping (ε, θ) 7→ Rε (θ, P b N ) = sup ℓ Aθ (z) , lim Rε (θ, P ε→∞
z∈Z
b N ) is monotonically nondecreasing for each θ ∈ Θ. Then according to Dini’s and ε 7→ Rε (θ, P b N ) converges uniformly for ε → ∞, i.e., theorem, Rε (θ, P b lim sup sup ℓ Aθ (z) − Rε (θ, PN ) = 0. ε→∞ θ∈Θ
z∈Z
As Θ is compact, supθ∈Θ supz∈Z ℓ(Aθ (z)) exists and is an ε-independent term. Therefore, uniform convergence allows interchanging the limit and the infimum: b N ) = lim inf Rε (θ, P b N ) = inf lim Rε (θ, P b N ) = inf sup ℓ Aθ (z) , lim Rε (θε⋆ , P ε→∞
ε→∞ θ∈Θ
θ∈Θ ε→∞
θ∈Θ z∈Z
which completes the proof.
B.3
■
Proof of Theorem 2
We now prove Theorem 2. Proof of Theorem 2. According to Lemma 6, b N ) ≤ Rε (θ, P b N ) ≤ sup ℓ Aθ (z) , R(θ, P z∈Z
b N ), we get holds true for all θ ∈ Θ. From θε⋆ ∈ argminθ∈Θ Rε (θ, P b N ) ≤ Rε (θ⋆ , P b N ) ≤ sup ℓ Aθ⋆ (z) = inf sup ℓ Aθ (z) . Rε (θε⋆ , P PEP PEP θ∈Θ z∈Z
z∈Z
Taking inf θ∈Θ consecutively from smaller to larger terms, we get b N ) ≤ inf Rε (θ, P b N ) = Rε (θ⋆ , P b N ) ≤ inf sup ℓ Aθ (z) . inf R(θ, P ε
θ∈Θ
θ∈Θ z∈Z
θ∈Θ
Applying Proposition 3 with Y = Aθ (z), bN) = Rε (θ, P
sup
E (ℓ(Aθ (z))) ≤ E
b N ) z∼Q Q∈Uε (P
bN z∼P
32
ℓ(Aθ (z)) + ε Lip(ℓ),
holds for every θ ∈ Θ. Take inf θ∈Θ on the left-hand side of the inequality to have b N ) = inf Rε (θ, P b N ) ≤ E ℓ(Aθ (z)) + ε Lip(ℓ), Rε (θε⋆ , P θ∈Θ
bN z∼P
and also take inf θ∈Θ on the right-hand side to get b N ) ≤ inf E Rε (θε⋆ , P
θ∈Θ z∼P bN
ℓ(Aθ (z)) + ε Lip(ℓ).
b N ) with We conclude the proof by upper-bounding true risk R(θε⋆ , P) by DRO risk Rε (θε⋆ , P probability at least 1 − β, using Theorem 1. ■
33
C
Additional details of Section 4.1
bi , and Fbi are conRemark 3 (Well-posedness of (DR-L2O)). The entries of Am (θ), bm (θ), G tinuous in θ, so the data of (DRO-PEP-D) vary continuously with θ. By Lemma 2, the Slater condition holds for (DRO-PEP-D) for every ε > 0, and compactness of Θ (Assumption 2) b N ) is continuous by Berge’s maximum implies it holds uniformly over Θ. Hence θ 7→ Rε (θ, P theorem [10, Chapter 6], and problem (DR-L2O) attains its minimum over Θ. Algorithm 1 describes the solution method for the learning problem (DR-L2O). Algorithm 1 Stochastic Gradient Method for Distributionally-robust L2O (DR-L2O) bN , parameter initialization θ0 , mini-batch size n, iteration Input: training data D counter k = 0. repeat bn from training data D bN Sample mini-batch D b n of D bn Define mini-batched empirical distribution P b Evaluate dRε (θ, Pn )/dθ using Lemma 7. Update θk+1 ← AdamWUpdate(θk ). k ← k + 1. until θk converges b n ) given a dataset D bn of n samples, we use In order to evaluate the gradient of θ 7→ Rε (θ, P the approach of [2, 1, 32] to differentiate the solution mapping of conic linear problem in terms b of problem parameters. Let T S, Dn be a solution mapping of the problem (DRO-PEP-D) b N ) is and let L be an objective function of (DRO-PEP-D). Then the DRO risk Rε (θ, P composed of solution mapping and L as follows: b n ) = L T Sθ , Aθ (D bn ) . Rε (θ, P We finally apply the chain rule to obtain the gradient with respect to θ, as the following: b n ) evaluates as bn , the gradient of θ 7→ Rε (θ, P Lemma 7. Given D d bn ) = dL θ 7→ Rε (θ, D dθ dt t=T (Sθ ,Aθ (Dbn )) ! M n X ∂T d(Am (θ), bm (θ)) X ∂T dAθ (ẑi ) + . × bi , Fbi ) ∂(A , b ) dθ dθ b m m ∂( G (S ,A ( D )) n θ θ b m=0 i=1 (Sθ ,Aθ (Dn ))
b n ) as the parametric conic program (DRO-PEP-D) with Proof. We interpret the risk Rε (θ, P b b parameters (Gi , Fi ) for i = 1, . . . , n and (Am , bm ) for m = 0, 1, . . . , M . Using the chain rule on b n ) = L T Sθ , Aθ (D bn ) , Rε (θ, P gives the desired result.
■ 34
1.0
wk
0.8 0.6 0.4 0.2
1
5
10
15
k Figure 6: A plot showing the exponential weight decay for our training objective with weights wk = 0.9K−k for K = 15.
Training procedure. For each experiment, we fix the algorithm structure and learn the step-size schedule θ = {θk }K−1 k=0 for each horizon K and each of (DR-L2O), (OPT-PEP), and (L2O). For (DR-L2O) and (L2O), we use minibatched stochastic gradient descent with the AdamW optimizer [40, 50], with weight decay cross-validated over {0, 10−5 , 10−4 , 10−3 }, combined with a linear warm-up followed by a cosine annealing learning rate schedule [49]. For (OPT-PEP), we follow the approach of [37]: their [37, Theorem 3.2] gives an explicit gradient of the worst-case rate with respect to the step size, and we apply gradient descent with the same cosine annealing schedule. All experiments run for 1000 iterations, with 10% of the iterations used for a linear warm-up to the maximum learning rate, which is cross-validated from {10−5 , 10−4 , 10−3 }. When learning the schedule, we do not use the performance loss at the last iterate as the training objective as this can hinder gradient information flow to earlier stepsizes [18]. Instead, we use an exponentially decaying weighted summation for the training loss. Letting ℓk represent the loss at iterate k, we minimize the sample analogue of the surrogate objective: "K # X E wk ℓk (Aθ (z)) . z∼P
k=1
We choose wk = 0.9K−k , and show a plot of the exponential decay in Figure 6. Additionally, during the learning procedure, in order to remove numerical issues with any step size component becoming negative, we square the step size before plugging it into the first order method. We initialize the values so that the square of the input step size is our desired value but we backpropagate the gradients to the square root. Parameter selection. For each experiment, we set a maximum value of K and learn a separate schedule for each k = 1, . . . , K used as the learning horizon. The Wasserstein 35
radius ε > 0 from Section 3.2 controls the probability of constraint satisfaction; in practice it is treated as a hyperparameter tuned by cross-validation [29, 61]. We cross-validate over ε ∈ {10−2 , 10−1 , 1, 5, 10} on a validation set and choose the schedule minimizing the empirical risk. A separate test set is used for final loss evaluations.
36
D
Further experiment details
All examples are written in Python 3.13 with JAX version 0.9.0; the code is available at https://github.com/stellatogrp/dro_pep. The SDPs are solved using the Clarabel solver [30] with primal residual, dual residual, and duality gap tolerances all set to 10−5 . For the derivatives, we use the diffcp package [1, 2] in order to differentiate through the KKT conditions of the inner SDP. All computations were run on a high performance computing cluster with 4 CPU cores and 20GB of RAM. In all plots, shaded bars show the 10th to 90th quantile ranges. Furthermore, in all SGD iteration timing tables, our code with JAX uses just-in-time (JIT) techniques for computational efficiency, and we noticed that the first few SGD iterations take significantly longer because of the compilation time. Since this is an initial startup cost, these initial times become outliers for the full distribution of times so we discard the first 5 iteration times before computing the averages and standard deviations. In the LASSO and TV inpainting experiments, since the optimal values are not always at 0 with objective value 0, we use a separate set of 1000 instances, solve all problem instances to optimality, and compute the maximum distance to optimality from the initial points. We then add a 10% buffer to ensure that all sample values in the training-set satisfy the interpolation conditions with this computed sample distance.
D.1
Further details of unconstrained quadratic experiment in Section 6.1
For each sample index i, we sample Qi from the Marčenko-Pastur distribution [53, 62] and zi0 from a uniform distribution over a ball of radius R around the origin. We rejection sample Qi ’s with λ(Qi ) ̸⊆ [µ, L]. For this experiment, the in-distribution parameters are d = 300, µ = 1, L = 10, R = 10. For the out-of-distribution set, we increase to L = 11 and hold all other parameters constant. We use a training set of size 1000; the validation set and both in-distribution and outof-distribution test sets are size 250. We use a mini-batch of size 20 for the SGD iterations. For the initialization, we initialize the stepsize schedules with constant value 1.5/(µ + L). In Figure 7 we show the test losses of all the learned schedules while in Table 1 we provide average times per SGD iteration for select values of K.
D.2
Further details of LASSO experiment in Section 6.2
We use a sparse coding example, or the problem of recovering a sparse vector x̃ from noisy measurements through a sparse dictionary matrix A ∈ Rm×n [20]. We form A by sampling each entry Aij ∼ N (0, 1/m) and normalizing each column to have unit 2-norm; we fix A for both the in-distribution and out-of-distribution sets and compute L as the largest eigenvalue of AT A. We sample a sparse vector x̃j ∈ Rn as x̃j ∼ N (0, σx2 ) with probability pmask and 0 2 otherwise. We then construct bj = Ax̃j + εj where εj ∼ N (0, σerr ). In the experiment, we 37
In-distribution
Out-of-distribution
Avg. f (xK ) − f (x? )
101 100
100 10−1
10−2
10−2 10−3
10−4
10−4 2
4
6
8
10
12
14
2
4
6
K
8
10
12
14
K
L2O
DR-L2O
OPT-PEP
Figure 7: Quadratic minimization experiment (Section 6.1) test losses across horizon K on a Left: in-distribution test set and Right: out-of-distribution test set.
set (m, n) = (250, 500), λ = 0.4, σx = 2.0 for the in-distribution set and σx = 3.0 for the out-of-distribution set, σerr = 0.01, and pmask = 0.1. Similarly to the unconstrained quadratic experiment, the training set, validation set, and test set sizes are 1000, 250, and 250, respectively. We use a size 10 mini-batch for SGD and initialize the stepsize schedules with constant value 1/L. We initialize all x0 values at 0. In Figure 8 we show the test losses of all the learned schedules while in Table 2 we provide average times per SGD iteration for select values of K.
D.3
Further details of TV inpainting experiment in Section 6.3
For the in-distribution dataset, we use the Olivetti faces dataset [69, 14], which contains 40 faces and 10 grayscale images each (64 × 64 pixels) of different angles. For the corruption, we choose 10% of pixels of each image uniformly at random to black out (i.e. set to 0). We randomly pick 28 faces for the training set, 4 faces for the validation set, 8 faces for the test set, and ensure that the images are separated in a stratified manner to prevent data leakage during training. For the out-of-distribution set, we pick 40 random color images from Tiny ImageNet [43, 22], where each Ui,j ∈ R3 is now a triplet of RGB values. We reformulate the TV LP in the standard LP form minimize cT x subject to Ax = b Gx ≤ h x ≤ x ≤ x, matching the formulation used in [52, Equation 1] and [11], and its saddle-point reformu38
Table 1: Quadratic minimization experiment per-SGD-iteration wall time (mean ± 2σ, in seconds) for the validation-best schedule used by the losses figure. Framework
K
Time (s)
L2O
1 5 10 15
0.003 ± 0.001 0.006 ± 0.001 0.009 ± 0.005 0.011 ± 0.001
DR-L2O
1 5 10 15
1.181 ± 1.102 1.740 ± 1.135 1.969 ± 0.360 11.064 ± 1.261
OPT-PEP
1 5 10 15
0.009 ± 0.004 0.017 ± 0.006 0.050 ± 0.005 0.202 ± 0.011
lation. To obtain this form from the TV LP in Section 6.3, we vectorize the image U into vec(U ) ∈ Rmn and introduce nonnegative slacks t ∈ R2(m−1)(n−1) that upper-bound each absolute-difference term in the objective. Let D stack the horizontal and vertical pixeldifference operators and E select the known-pixel coordinates. Writing x = (vec(U ), t), the data of the standard LP are " # " # h i 0 D −I c= , A = E 0 , b = vec(U orig )K , G = , h = 0, 1 −D −I with box bounds x = 0 and x = (1, +∞) on (vec(U ), t). The corresponding PDHG iterations are xk+1 = proxτ k f xk − τ k M T uk , x̄k+1 = xk+1 + ρk (xk+1 − xk ),
uk+1 = proxσk g∗ uk + σ k M x̄k+1 . T
h
T
T
i
For each sample LP, we define the stacked constraint matrix M = A −G and compute the maximum value of ∥M ∥2 across the samples, denoted Mmax . We initialize the training procedures with (τ k , ρk , σ k ) = (0.5/Mmax , 1, 0.5/Mmax ). In order for the Lagrangian duality gap to make sense as a metric, we need to initialize the algorithm from a feasible point, so we choose to initialize x0 = 1/2 and u0 = 1. In Figure 9 we show the test losses of all the learned schedules and in Figure 10 we show the fractions of problems solved. In Figure Lastly, in Table 3 we provide average times per SGD iteration for select values of K.
39
In-distribution
Out-of-distribution
103
Avg. f (xK ) − f (x? )
101 102 100
101
10−1
100
10−2
10−1 2
4
6
8
10
12
14
2
4
6
K L2O-ISTA
8
10
12
14
K L2O-ALISTA
DR-L2O
OPT-PEP
Figure 8: LASSO experiment (Section 6.2) text losses against horizon K on a Left: in-distribution test set and Right: out-of-distribution test set.
In-distribution
Out-of-distribution
Avg. Lagrangian gap
102
102 101
2
4
6
8
10
2
4
K
6
8
10
K L2O
DR-L2O
OPT-PEP
Figure 9: PDLP experiment (Section 6.3) test losses against horizon K on the Left: Olivetti in-distribution test set and Right: Tiny ImageNet out-of-distribution test set.
40
Test set, fraction of problems solved
Out-of-distribution
In-distribution
η = 0.01
η = 0.05
η = 0.1
1
1
1
0
0
0
1
1
1
0
0 2
4
6 K
8
10
L2O
0 2
4
6 K
8
10
DR-L2O
2
4
6 K
8
10
OPT-PEP
Figure 10: Fractions of PDLP problems solved to different relative tolerances against horizon K. Top: Fractions solved on the Olivetti in-distribution set. Bottom: Fractions solved on the Tiny ImageNet out-of-distribution set.
Original
Corrupted
Opt. LP Reconstruction
OPT-PEP
L2O
DR-L2O
Figure 11: More reconstructions on images from the Tiny ImageNet dataset. The images are ordered in the same way as in Figure 5.
41
Table 2: LASSO experiment per-SGD-iteration wall time (mean ± 2σ, in seconds) for the validation-best schedule used by the losses figure. Framework
K
Time (s)
L2O-ISTA
1 5 10 15
0.006 ± 0.000 0.015 ± 0.001 0.025 ± 0.002 0.037 ± 0.002
L2O-ALISTA
1 5 10 15
0.004 ± 0.000 0.011 ± 0.001 0.020 ± 0.001 0.029 ± 0.001
DR-L2O
1 5 10 15
0.144 ± 0.027 0.825 ± 0.129 2.993 ± 0.327 14.697 ± 8.933
OPT-PEP
1 5 10 15
0.015 ± 0.001 0.050 ± 0.009 0.364 ± 0.041 2.008 ± 0.424
Table 3: TV inpainting experiment per-SGD-iteration wall time (mean ± 2σ, in seconds) for the validation-best schedule used by the losses figure. Framework
K
Time (s)
L2O
1 5 10
0.119 ± 0.014 0.450 ± 0.016 0.845 ± 0.019
DR-L2O
1 5 10
0.410 ± 0.018 3.755 ± 0.154 57.125 ± 27.471
OPT-PEP
1 5 10
0.037 ± 0.012 0.371 ± 0.021 3.954 ± 0.309
42