Learning to Relax: Setting Solver Parameters Across a Sequence of Linear System Instances
Authors: Mikhail Khodak
Edmond Chow
Maria Florina Balcan
Ameet Talwalkar
Venue: ICLR 2024 spotlight (2024)
Keywords: scientific computing
data-driven algorithm design
online learning
multi-armed bandits
contextual bandits
numerical analysis
learning-augmented algorithms
algorithms with predictions
Primary area: learning theory
TL;DR: We show provable guarantees for learning the relaxation parameter of linear system solvers.

Abstract: Solving a linear system ${\bf Ax}={\bf b}$ is a fundamental scientific computing primitive for which numerous solvers and preconditioners have been developed. 
	These come with parameters whose optimal values depend on the system being solved and are often impossible or too expensive to identify;
	thus in practice sub-optimal heuristics are used.
	We consider the common setting in which many related linear systems need to be solved, e.g. during a single numerical simulation.
	In this scenario, can we sequentially choose parameters that attain a near-optimal overall number of iterations, without extra matrix computations?
	We answer in the affirmative for Successive Over-Relaxation (SOR), a standard solver whose parameter $\omega$ has a strong impact on its runtime.
	For this method, we prove that a bandit online learning algorithm—using only the number of iterations as feedback—can select parameters for a sequence of instances such that the overall cost approaches that of the best fixed $\omega$ as the sequence length increases.
	Furthermore, when given additional structural information, we show that a _contextual_ bandit method asymptotically achieves the performance of the _instance-optimal_ policy, which selects the best $\omega$ for each instance.
	Our work provides the first learning-theoretic treatment of high-precision linear system solvers and the first end-to-end guarantees for data-driven scientific computing, demonstrating theoretically the potential to speed up numerical methods using well-understood learning algorithms.

LEARNING TO RELAX: SETTING SOLVER PARAMETERS ACROSS A SEQUENCE OF LINEAR SYSTEM INSTANCES
Mikhail Khodak
Edmond Chow
Maria-Florina Balcan
Ameet Talwalkar
CMU
Georgia Tech.
CMU
CMU
khodakk cmu edu
cchow@cc.gateche .dd
ninamff C8s cmu edu
Galwalkar@cmm edh.
ABSTRACT
Solving a linear system Ax b is a fundamental scientific computing primitive for which numerous solvers and preconditioners have been developed. These come with parameters whose optimal values depend the system being solved and are often impossible or too expensive to identify; thus in practice sub-optimal heuristics are used. We consider the common setting in which many related linear systems need to be solved, e.g. during a single numerical simulation. In this scenario, can we sequentially choose parameters that attain a near-optimal overall number of iterations, without extra matrix computations? We answer in the affirmative for Successive Over-Relaxation (SOR), standard solver whose parameter has a strong impact on its runtime. For this method, we prove that a bandit online only the number of iterations as feedback can select parameters for sequence of instances such that the overall cost approaches that of the best fixed w as the sequence length increases. Furthermore, when given additional structural information, we show that contextual bandit method asymptotically achieves the performance of the instance-optimal policy, which selects the best w for each instance. Our work provides the first carning-theorette treatment of sish-precision linear system solvers and the first end-to-end guarantees for data-driven scientific computing, demonstrating theoretically the potential to speed up numerical methods using well-understood learning algorithms.
learning algorithm-
INTRODUCTION
The bottleneck subroutine in many scientific computations is a solver returning approximate so- lution to linear system. For example, simulating a partial differential equation (PDE) often involves solving sequences of nigh-dimensional systems to very high precision (Thomas, 1999). A vast array of solvers and preconditioners have thus been developed, many of which have tunable parameters that significantly affect runtime (Greenbaum, 1997; Hackbusch, 2016). There is long literature analyzing these algorithms, and indeed for some problems we have strong understanding of the optimal parame- ters for a given matrix. However, computing them can be more costly than solving the original system, leading to an assortment of heuristics for setting good parameters (Ehrlich, 1981; Golub & Ye, 1999).
We provide an alternative to such heuristics by taking advantage of the fact that we often sequentially solve many linear systems. In addition to numerical simulation, this occurs in graphics computations such as mean-curvature flow (Kazhdan et al., 2012), nonlinear system solvers (Marquardt, 1963), and beyond. natural approach is to treat these instances as data to passed to a machine learning (ML) algorithm; in particular the framework of online learning (Cesa- (Cesa-Bianchi & Lugosi, 2006) provides language to reason about such sequential learning problems. For example, we otherwise would solve sequence of linear systems (A1, bj)...., (AT, br) using given solver with fixed parameter, can we use ML to do as well as the best choice of that parameter, i.e. can we minimize regrer? Or, if the matrices are all diagonal shifts of single matrix A, can we learn the functional relationship between the shift e and the optimal solver parameter for A, = A + c.Ip. i.c. can we predict using context?
We investigate these questions for the Successive Over-Relaxation (SOR) solver, a generalization of Gauss-Seidel whose relaxation parameter Figure 1, noting the log-scale). SOR and its symmetric variant are well-studied and often used as pre- conditioners for Krylov methods such as conjugate gradient (CG), bases for gemi-iterative schemes, and as multigrid smoothers. Analogous to some past setups in data-driven algorithms (Balcan et al., 2018; Khodak et al., 2022), we sequentially set the parameter w, for SOR to use when solving each lin-
(0, 2) dramatically affects the number of iterations (c.f.

---

system (A,, b.). Unlike past theoretical studies of related methods (Gupta & Roughgarden, 2017;
Bartlett et al., 2022; Balcan et al., 2022), we aim to provide end-to-ens end guarantees covering the full pipeline from data-intake to efficient learning to orccution-while minimizing dependence on the di- mension (n can be 10 or higher) and precision (1 /s can be 10 or higher). We emphasize that we do not seek to immediately improve the empirical state of the art, and also that existing research on saving computation when solving sequences of linear systems (recycling Krylov subspaces, reusing precon- ditioners, etc.) complementary to our own, i.c. can be used in addition to the ideas presented here.
1.1 CORE CONTRIBUTIONS
We study two distinct theoretical settings, corresponding to views the problem from two different approaches to data-driven algorithms. In the first we have a deterministic sequence of instances and study the spectral radius of the iteration matrix, the main quantity of interest in classical analysis of SOR (Young, 1971). We show how to convert its asymptotic guarantee into a surrogate loss that upper bounds the number of iterations via a quality measure of the chosen parameter, in the style of algorithms with predictions (Mitzenmacher & Vassilvitskii, 2021). The bound holds under near- asymptotic condition implying that convergence occurs near the asymptotic regime, i.e. when the spec- tral radius of the iteration matrix governs the convergence. We verify the assumption and show that one can learn the surrogate losses using only bandit feedback from the original costs; notably, despite being non-Lipschitz, we take advantage of the losses' unimodal structure to match the optimal ((T2/3) re- gret for L Lipschitz bandits (K (Kleinberg, 2004). Our bound also depends only logarithmically on the preci- sion and not at all on the dimension. Furthermore, we extend to the diagonally shifted setting described before, showing that an efficient, albeit pessimistic, contextual bandit (CB) method has ((T3/4) re- gret w.r.t. the instance-optimg policy that always picks the best w/. Finally, we show similar analysis of learning a relaxation parameter for the more popular (symmetric SOR cOR-precooniitione) CG method.
Our second setting semi- stcchastchas,ic, with target vectors b, drawn i.i.d. from a (radially truncated) Gaussian. This is reasonable simplification, as convergence usually depends more strongly on Ar, on which we make no extra assumptions. We show that the expected cost of running a symmetric variant of SOR (SSOR) is optimal number of iterations-raaie rather than with the best upper bound-ame (b) analyze more practical, regression- -ased CB algorithms (Foster & Rakhlin, 2020; Simchi-Levi & Xu, 2021). We then show Ö (/TTa/n) regret when comparing to the single best w and (TT9/11 Vn) regret w.r.t. the instance-optimal policy in the diagonally shifted setting using a novel, Chebyshev regression-basee CB algorithm. While the results do depend on the dimension n, the dependence is much weaker than that of past work on data- tuning of a related regression problem (Balcan et al., 2022).
2 )Lipschitz w.r.t. w, so we can (a) compete with the
Remark 1.1. Likely the most popular algorithms for linear systems are Krylov subspace methods such as CG. While an eventual aim of our line of work is to understand how to tune (many) parameters of (preconditioned) CG and other algorithums, SOR is as a meaningful starting point. In fact, we show that our near-aste eear-asymptotic analysis extends directly, and in the cami--tochastic setting there is natural path to (e.g.) SSOR- preconditioned CG, as SSOR-preconditioner can be viewed as computing polynomials of iteration matrices where SSOR just takes powers. Lastly, apart from its use as preconditioner and smoother, SOR is still sometimes preferred for direct use as well (Fried & Metzler, 1978; Van Vleck & Dwyer, 1985; King et al., 1987; Wožnicki, 1993; 2001).
method and serves
1.2 TECHNICAL AND THEORETICAL CONTRIBUTIONS
By studying a scientific computing problem through the lens of data-driven algorithms and online learning, we also make the following contributions to the latter two fields:
1. Ours is the first head-to-head comparison of two leading theoretical approaches to data-driven algo- rithms applied to the same problem. While the algorithms with predictions approach in Section takes better advantage of the scientific computing literature to obtain (arguably) more interpretable and caneension-indepenee with the quantity of interest in Section 3 and enables guarantees for modern CB algorithms.
bounds, data-driven algorithm design (Balcan, 2021) competes directly
2. For algorithms with predictions, our near-asymptotic approach may be extendable to other iterative solvers, we demonstrate with CG. We also show that such performance bounds (partially- observable) cost are learnable even when the bounds themselves are too expensive to compute. 3. In futa-driven algorithm design, we take the novel theoretical approach of proving continuity of the expectation of a discrete cost, rather than showing dispersion of its discontinuities (Balcan et al., 2018) or bounding predicate complexity (Bartlett et al., 2022).
2

---

Algorithm 1: Successive over-relaxation (SOR) with relative convergence condition.
Input: A E Rax", b E IR", parameter E 2), initial vector E R", tolerance > 0 // D is diagonal, is strictly lower triangular // compute the third normal form // compute initial residual
D + L 1++D ++L* + A
W D/u
To < Ax
for k 0.... do
|| | S <e||rolla then
L return K
// return iteration count (for use in learning) // update vector after solving triangular system // compute the next residual
= W Ik
Tk+1 - b Ax
4.
We introduce the idea of using CB to set mnstance-adaptii algorithmic parameters; while (linear) instance-adaptivity was also shown via convexity by Khodak et (2022), we go further by taking advantage of multi-instance structure to asymptotically do as well as the instunce-optinal policy. We show that standard discretization-bant -based bandit algorithms are optimal for sequences of adver- sarially chosen semi- semi-Lipschitz losses that generalize regular Lipschitz functions (c.f. Appendix B). We introduce a new CB method that combines SquareCB (Foster & Rakhlin, 2020) with Chebyshev polynomial regression to get sublinear regret Lipschitz losses (c.f. Appendix C).
1.3 RELATED WORK AND COMPARISONS
We discuss the existing literature solving sequences of linear systems (Parks et al., 2006; Tebbens & Tûma, 2007; Elbouyahyaoui et al., 2021), work integrating ML with scientific computing to amortize cost (Amos, 2023; Arisaka & Li, 2023), and past theoretical studies of data-driven algorithms (Gupta & Roughgarden, 2017: Balcan et al., 2022) in Appendix A. For the latter we include a detailed com- parison of the generalization implications of our work with the GJ framework (Bartlett et al., 2022). Lastly, we address the baseline of approximating the spectral radius of the Jacobi iteration matrix.
2 ASYMPTOTIC ANALYSIS OF LEARNING THE RELAXATION PARAMETER
We start this section by going over the problem setup and the SOR solver. Then we consider the asymptotic analysis of the method to derive a reasonable performance upper bound to target as a surrogate loss for the true cost function. Finally, we prove and analyze online learning guarantees.
2.1 SETUP
At each step 1 1,..., T of (say) numerical simulation we get a linear system instance, defined by matrix-vector pair (A,, b,) E Raxa R., and asked for a vector x E R such that the norm of its residual or defect b, - A is small. For now we define "small" in relative sense, specifically |A,x-b;lla byll2 S =||b:]l2 for some tolerance e E (0, 1); note that when using an iterative method initial- ized at x = 0, this corresponds to reducing the residual by factor 1/e, which we call the precision. In applications can be quite high, and so we will show results whose dependence on is at worst logarithmic. To make the analysis tractable, we make two assumptions (for now) about the matrices A: they are symmetric positive-definite and ansistently-ordeees tion 4.23)). We emphasize that, while not necessary for convergence, both are standard in the analysis of SOR (Young, 1971); see Hackbusch (2016, Criterion 4.24) for multiple settings where they holds.
......,T
(c.f. Hackbusch (2016, Defini-
To find suitable x for each instance in the sequence we apply Algorithm 1 (SOR), which a high-level works by multiplying the current residual r by the inverse of matrix der-deiveed from the diagonal D and lower-triangular component L of A and then adding the result to the current iterate x. Note that multiplication by W_ is efficient because W., is triangular. We will measure the cost of this algorithm by the number of iterations takes to reach convergence, which we denote by SOR(A. A, b, w), or SOR((W) for short when it is run on the instance (A,, b.). For simplicity, we will assume that the algorithm is always initialized at x = O,,, and so the first residual is just
Having specified the computational setting, we now turn the learning objective, which to sequentially set the parameters Wi...., ......,,,y w. so as to minimize the total number of iterations:
SOR((wr) = Σ' t=1
SOR(A,, wr)
(1)

---

To set Wt at some time t> 1, we allow the learning algorithm access to the costs (OR,(w,) incurred at the previous steps 8 = 1....,t - 1: in the literature on online learning this is referred to as the bandit or partial feedback setting, to distinguish from the (easier, but unreasonable for us) full information case where we have access to the cost function SCR, at every w in its domain.
Selecting the optimal w; using no information about A, is impossible, so we must use a comparator to obtain an achievable measure of performance. In online learning this is done by comparing the total cost incurred (1) to the counterfactual cost had we used a single, best-in- -iinddsihhiihht w at every timestep t. We take the minimum over some domain 8 C (0, 2). as SOR diverges outside it. While in some settings we will compete with every w E (0, 2). we will often algorithmically use [1, Waax] for some Wmax < 2. The upper limit ensures bound on the number of iterations required by bandit algorithms-ang and the lower limit excludes w < 1, which is rarely used because theoretical convergence of vanilla SOR is worse there for realistic problems, e.g. those satisfying our assumptions.
This comparison-bar based approach for measuring performance is standard in online learning and effectively assumes good w E Ω that does well-enough on all problems; in Figure (center-left) we show that this is sometimes the case. However, the center-right plot in the same figure shows might do better by using additional knowledge about the instance; in online learning this is termed a context and there has been extensive development of contextual bandit algorithms that do well the best fixed policy mapping contexts to predictions. We will study an example of this in the diagonally shifted setting, which A, = A+I,,, for scalars C E I; while mathematically simple, this structure arises in natural settings, e.g. solving the heat equation with temporally variable diffusivity, and is well-motivatee by other applications (Frommer Glässner, 1998; Bellavia et al., 2011; Baumann & van Gijzen, 2015; Anzt et al., 2016; Wang et al., 2019). Furthermore, the same learning algorithms can also be extended to make use of other context information, e.g. rough spectral estimates.
2.2 ESTABLISHING SURROGATE UPPER BOUND
Our first goal is to solve T linear systems almost as fast if we had used the best fixed w Ω. In on- line learning, this corresponds to minimizing regret, which for cost functions & : S2 + R is defined as
(egret,2(44fFe = =E (e(wer) min
21=1 (a(w) WESt
T
(2)
In particular, since we can upper-bound the objective (1) by cost min Et= (w). if we show that regret is sublinear in T then the leading-order term in the upper bound corresponds to the cost incurred by the optimal fixed w.
plus the optimal
Many algorithms attaining sublinear regret under different conditions on the losses E, have been developed (Cesa-Bianchi & Lugosi, 2006; Bubeck & Cesa-Bianchi, 2012). However, few handle losses with discontinuities-ie i.e. most algorithmic costs and and those that do (necessarily) need additional conditions on their locations (Balcan et al., 2018; 2020). At the same time, numerical analysis often deals more directly with continuous asymptotic surrogates for cost, such as convergence rates. Taking inspiration from this, and from the algorithms with predictions idea of deriving surrogate loss functions for algorithmic costs (Khodak et al., 2022), in this section we instead focus on finding upper bounds U, on SOR, that are both (a) learnable and (b) reasonably tight in in-practice. We can then aim for overall performance nearly as good as the optimal w E Ω as measured by these upper bounds:
T
(wr) S Σ Ui(w)) Regret ((U,)T + min wEll
Σ -u((w)=((T o(T) + min Ue(w) (3)
WEB 1=1
t=1
A natural approach get a bound U, is via the defect reduction matrix C. I, so named because the residual at iteration k is equal to C b and b is the first residual. Under our assumptions on A, Young (1971) shows that the spectral radius p(Cw) of C. is (nontrivial com- pute) piecewise function of w with a unique minimum in [1, 2). Since we have error b|/ |lbl2 S limn--oc VI[C 2 asymptotically bounds how much the error is re-
AD
||C 2l at iteration k. p(C..) duced at each step. It is thus often called the asymptotic convergence rate and the number of iterations is said to be roughly bounded by = Tog ((C.., (e.g. Hackbusch (2016, Equation 2b3. However, while it is tempting to use this as our upper bound U, in fact it may not upper bound the number of iterations at all. since C.,, is not normal and so in practacee the iteration often goes through a transient phase where the residual norm first increases before decreasing (Trefethen & Embrec, 2005, Figure 25.6).
log
Thus we must either take a different approach or make some assumptions. Note that one can in-fact show an w-dependent, finite-time convergence bound for SOR via the energy norm (Hackbusch,
4

---

amall ond bound rest-sttg metiii louund
MoSi vay bouadd
BRIItOr 2
0.0
Inattre
032
1.2
14 1.5 1.
12
12 18
Figure 1: Left: comparison of different cost estimates. Center-left: mean performance of different parameters across forty instances of form A+ 12c-31 In. where ~ Beta(2, (2 6). Center-right: the same but for e ~ Beta(1 1/2, which is relatively higher-variance. In both cases the dashed line indicates instance-optimal performance, the matrix A is discrete Laplacian of a 100 x 100 square domain, and the targets are truncated Gaussians. Right: asymptocity as measured by the difference between the spectral norm at iteration & and the spectral radius, together with its upper bound r(1 ((Cw))
2016, Corollary 3.45), but this can give rather loose upper bounds on the number of iterations (c.f. Figure (left). Instead, we make the following assumption, which roughly states that convergence always occurs near the asymptotic regime, where nearness is measured by parameter (0, 1):
Assumption 2.1. There exists E 10 s.I. A w E Ωt the matrix C. = I, ADD/w L)-1 satisfies ||C* ll S (p(C.) + 7(1 ()C
k = min | bll2"
bel2<ce|
This effectively assumes an upper bound p (C.... + r(1 ((Cw)) on the empirically observed convergence rate, which gives us a measure of the quality of each parameter w for the given instance (A, b). Note that the specific form of the surrogate convergence rate was chosen both because is convenient mathematically- is convex combination of I and the asymptotic rate pl P(Cw) ((C))-ade and because empirically we found the degree of "asymptocity" as measured by ||C* ||/ p(C..) for right before convergence vary reasonably similarly a fraction of - p( CC,.) (c.f. Figure (right)). This makes intuitive sense, as the parameters w for which convergence is fastest have the least time to reach the asymptotic regime. Finally, note that since limg- ||C*1l// = p( for every there always exists k' s.t. C l2 < (p(C..) + 7) A k 2 K'; therefore, since - p(C,) > 0, we view Assumption 2.1 not as a restriction on C. (and thus on A). but rather as an an assumption on e and b. Specifically, the former should be small enough that C reaches that asymptotic regime for some before the criterion bll2 I e|bbll is met; for similar reasons, the latter should not happen to be an eigenvector corresponding to tiny eigenvalue of C, (c.f. Figure 2 (left)).
0
Having established this surrogate of the spectral radius, we can use it to obtain a reasonably tight upper bound U on the cost (c.f. Figure (left)). Crucially for learning, we can also establish the following properties via the functional form of p(C..) derived by Young (1971):
Lemma 2.1. Define (ww =
logg (T(1))(())))))))) T (1 7) ma Wax 1}.
and w = /²/(1 /1 1/1-3211 32)2, where p(In D-1A). Then the following holds:
1. U bounds the number of iterations and itself bounded:: SOR(A A, b, w) < U(w) S - log
log
U is decreasing towards w", and -(1-r)logs log" -Lipschitz on w 2 w if 2 1 or 32 P 1
max Waax 1} that appears in the upper bounds on (ww) and in its Lipsehitz constant. This quantity will in some sense measure the difficulty of learning: if is close to for many of the instances under consideration then learning will be harder. Crucially, all quantities in the result are spectral and do not depend on the dimensionality of the matrix.
Lemma 2.1 introduces a quantity = T
2.3 PERFORMING AS WELL AS THE BEST FIXED W
Having shown these properties of U, we now show that it is learnable via Tsallis INF (Abernethy et al., 2015; Zimmert& & Seldin, 2021), a bandit algorithm which each instance L. samples we from discrete probability distribution over grid of dl relaxation parameters, runs SOR with w/ on the linear system (A,, by), and uses the number of iterations required SOR, (we) as feedback to update the probability distribution over the grid. The scheme is described in full in Algorithm 2. Note that it is a relative of the simpler and more familiar Exp3 algorithm (Auer et al., 2002), but has a slightly better dependence on the grid size d. In Theorem 2. we bound the cost of using the parameters W, suggested by Tsallis-INF INF by the total cost of using the best fixed parameter ws E Ω at all iterations as measured by the surrogate bounds U,-plus term that increases sublinearly in 7 and a term that decreases in the size of the grid.
5

---

Algorithm 2: Online tuning of a linear system solver using Tsallis-INE: The probabilities can be computed using Newton's method (e.g. Zimmert & Seldin (2021, Algorithm 2)).
Input: solver SOLVE: RRax" x R" x + Z,00 instance sequence {(A,, b.)), C RRRx xx""
CR
normalization K > 0, parameter grid g E 2 step-size 17 > 0
+ 0d for t 1, T do
// initialize vector of cumulative costs
p k arg minper (k, p)
// compute probabilities // sample action from grid ki.]] t kt.! + (SOLVE wr) - 1)/Pk) // run solver and update cost
4k Σ PP
sample in E [d] w.p. Pb. and set W/ = 8l4
Theorem 2.1. Define 04 + Te) max (B} Wime
1}, where 3e pI,, D A1) and Te
is the minimal T satisfying Assumption 2 and the second part of Lemma 2. If we run Algorithm using SOR initialized at x On as the solver, 8(i 1 + (Wmax 1) as the parameter grid,
normalization K 2 Tlog camax
Ing for O'max- = max; Ott, and step-size = VT then the expected
number of iterations is bounded as
T
T
logg
T Σ Us(w)
E
SOR((wr)
S 2K 2KT
+
min
(4)
d dlog ae
2
(0.
t=1
1=1
T Σ
Using Wann = 1+max
K
lloge log amax
and d I log" Cor 7
yields
t
1+/2--3 33
7
log"
1
27/T22
T Σ
3log 272
1 V
T Σ Ue(w) (5)
E
(w)) S 33log-
min WE (0.)
Ui(w)
+ min
8
$log*
amx
0mx
& Vlog's
4
O'max
(0,2)
Thus asymptotically (as T' oo) the average cost on each instance is that of the best fixed w E (0, 2). as measured by the surrogate loss functions Uz(w). The result clearly shows that the difficulty of the learning problem can be measured by how close the values of a are to one. As quantitative example, for the somewhat "easy" case of T. I 0.2 and B. I 0.9, the first term is < Tlog -.e. we take at most log 1L excess iterations average-aften around 73K instances.
The proof of Theorem 2.1 (c.f. Section E) takes advantage of the fact that the upper bounds U, are always decreasing wherever they are not locally Lipschitz: thus for any w E (0, Wans] the next highest grid value in g will either be better or (1/dd) worse. This allows us to obtain the same (T2333) rate as the optimal Lipschitz-bandi regret (Kleinberg, 2004), despite U, being only semi-Lipschitz. One important note is that setting Wime. K, and to obtain this rate involves knowing bounds on spectral properties of the instances. The optimal Wim requires a bound on max, 3r akin to that used by solvers like Chebyshev semi-iteration; assuming this and a reasonable sense of how many
iterations are typically required is enough to estimate Ommxx and then set d = V log" Orman" yielding
T/2 T/2_,
the right-hand bound in (5). Lastly, we note that Tsa adds quite little computational overhead: it has per-instance update cost of (dd) which for d O (T) is likely to be negligible in practice.
2.4 THE DIAGONALLY SHIFTED SETTING
The previous analysis is useful when a fixed w is good for most instances (A,, b.). A non-fixed comparator can have much stronger performance (c.f. the dashed lines in Figure 1 (center)), so in this section we study how to use additional, known structure in the form of diagonal shifts: at all E [T]] A; = A+ + I, for some fixed A and scalar ct. It is easy to see that selecting instance-dependent wt using the value of the shift is exactly the contextual bandit setting (Beygelzimer et al., 201 1), in which the comparator is a fixed policy f R. Ω that maps the given scalars to parameters for them. Here the regret is defined by Regret, y ((4))==) = Σin (wr) - Σ (f(c)). Notably, if f is the optimal mapping from C to w then sublinear regret implies doing nearly optimally at every instance. In our case, the policy w" minimizing U, is well-defined function of A, (c.f. Lemma 2.1) and thus of e, (Young, 1971); in fact, we can show that the policy is Lipschitz w.r.t. Ct (c.f. Lemma E. This allows us to use a very simple algorithm diintttmdsscretizingg the space of offsets into m intervals and run- ning Tsallis- -NF separately on each to obtain 33344) regret w.r.t. the instance-optimg policy w":

---

11d toa't
dell Fieely 2yr
Treblo 2011
Cwecd
12
85
rusiic
Figure Left: solver cost for b drawn from a truncated Gaussian v.s. b small eigenvector of C1.4- Center-left: cost to solve 5K diagonally shifted systems A, = A + 12c,-31 I, for C ~ Beta(2, 6). Center-right: total ASSOR-precoddiiinn
CG iterations taken while solving the 2D heat equation with time-varying diffusion coefficient (used context) on different grids, as function of the linear sys- tem dimension. Right: (smoothed) parameters chosen at each timestep of one such simulation, over- laid on a contour plot of the cost of solving the system at step / with parameter (c.f. Appendix G).
Theorem 2.2 (c.f. Theorem E.I). Suppose all offsets Ct lie in [Cmins Cmin + C] for some (Amin(D) mmin +1 for Ams as in Theorem 2.1. Amox Vi- Aviin +Di+
Cmin - min (A), and define I. = 1+3max
Then there is a discretization of this interval s.t. running Algorithm separately on each sequence of contexts in each bin with appropriate parameters results in expected cost
LAT 4 1
E 7 (wR)(wr) S V log Ommxx
2ACCLT
(6)
log dimux
t=1
Observe that, in addition to az, the difficulty of this learning problem also depends on the maximum spectral radius 9em. of the Jacobi matrices I,- D Ac(())((c)) via the Lipschitz constant L w*
2.5 TUNING PRECONDITIONES CONJUGATE GRADIENT
CG perhaps the most-used solver for positive definite systems; while it can be run without tuning, in practice significant acceleration can be realized via a good preconditioner such as (symmetric) SOR. The effect of w on CG performance can be somewhat distinct from that of regular SOR, requiring a separate analysis. We use the condition number analysis of Axelsson (1994, Theorem 7.17) to obtain an upper bound (w)) on the number of iterations required CG A, b, to solve a system. While the resulting bounds match the shape of the true performance less exactly than the SOR bounds (c.f. Figure 4), they still provide a somewhat reasonable surrogate. After showing that these functions also semi- semi-Lipsehitz (c.f. Lemma E.2), we can bound the cost of tuning CG using Tsallis-INF:
Theorem 2.3. Set pr =(D,A,"), Pmax = maxt Pt, /I vu= + ET= VHy and Kimas maxx R(A)). If min, b4 is a positive constant then for Algorithm using preconditioned CG as the solver there exists a parameter grid g E [2v2 + 2, Wear] and normalization K >> such that
CG((ur) = (
log
Amaax
T
VIL2
min WE(O,2)
U,(w)
log²
w/umax V/mmxx +1
Observe that the rate in T remains the same as for SOR, but the difficulty of learning now scales mainly with the spectral radii of the matrices D,A,
A STOCHASTIC ANALYSIS OF SYMMETRIC SOR
Assumption 2.1 in the previous section effectively encodes the idea that convergence will not be too quick for a typical target vector b. e.g. will not be a low-eigenvalue eigenvector of C.,, for some oth- erwise suboptimal w (e.g. Figure 2 (left)). Another way of staying in "typical regime" is randomness, which is what we assume in this section. Specifically, we assume that b, = m, u, V E [T]. where u, R" is uniform on the unit sphere and m² is x random variable with n degrees of freedom truncated to J0, n).. Since the standard n--iimensional Gaussian is exactly the case of untruncated m can be described as coming from a radially truncated normal distribution. Note also that the exact choice of truncation was done for convenience; any finite bound 2 n yields similar results.
We also make two other changes: (1) we study symmetric SOR (SSOR) and (2) we use an absolute con- vergence criterion, i.e. ||rx||2 E, not ||r||| el]ro||a. Symmetric SOR (c.f. Algorithm 8) is very
7

---

similar to the original, except the linear system being solved at every step is now symmetric: W.
2-w W.. 'W"- Note that the defect reduction matrix C., II, AW is still not normal, but it is non-orthogonally) similar to a symmetric matrix, A A//2 SSOR is twice expensive per-iteration, but often converges in fewer steps, and is commonly used base method because of its spectral properties (e.g. by the Chebyshev semi-iteration, c.f. Hackbusch (2016, Section 8.4.1).
3.1 REGULARITY OF THE EXPECTED COST FUNCTION
We can then show that the expected cost E,SSOR(A A, b, w) is Lipschitz w.n.t w (c.f. Corollary F.1). Our main idea is the observation that, whenever the error || C& b la falls below the tolerance E, random- ness should ensure that it does not fall close to the threshold that the error C bll2 of a nearby w' is not also below e. Although clearly related to dispersion (Balcan et al., 2018), here we study the behav- ior of a continuous function around a threshold, rather than the locations of the costs' discontinuities.
Our approach has two ingredients, the first being Lipschitzness of the error || [|C bibbll at each iteration k w.r.t. w, which ensures / bl12 E (e, O[ w/]] if ||C 2 S E < || " bll2. The second ingredient is anti-ooncentration, specifically that the probability that b2l22 lands in (e, O(w wDD) is O[w w'l). While intuitive, both steps are made difficult by powering: for high k the random variable || bll2 is highly concentrated because pl p(C,w) < 1; in fact its measure over the interval is O(w w'1 1lp((C... )*) To cancel this, the Lipschitz constant of |
w bil2 must scale with p p(C.)) which we can show because switching to SSOR makes Č is similar a normal matrix. The other algorithmic modification we make-using absolute rather than relative tolerance- is so that || w bll is (roughly) sum of i..d. x" random variables; note that the square of relative tolerance criterion | b [bbl does not admit such result. At the same time, absolute tolerance does not imply an a.s. bound on the number of iterations if ||b|22 is unbounded, which is why we truncate its distribution.
E, (SOR((())) can be bounded using Jensen's inequality by the probability that w and w' have different costs k o I, which is at most the probability that || bll2 or bll2 land in an interval of length O[K w'). Note that the Lipschitz bound in- cludes an O (yn) factor, which results from Č₁ having stable rank < due powering. Regularity of SSOR leads directly to regret guarantee for the same algorithm as before, Tsallis-INF:
Lipschitzness follows because [ESSSOR(w)
Theorem 3.1. Define Kmax max, (((A)) to be the largest condition number and 8min
min, pInn DD,' A.). Then there exists K 22log ") s.t. running Algorithmm with SSOR has regret
T
T
32K*T
2nn
E SSOR( (wi)
- min WE [
SSOR((W) S 2KV
84 miu d V
(8)
VTTVVn) Note that, while this shows convergence to the true optimal parameter, the constants in the regret term are much worse, not just due to the dependence be viewed as proof of the asymptotic (T + oo) correctness of Tsallis-NNF for tuning SSOR.
Setting d O(K2 VnTT yiclds regret bound of Ollog" 12
but also in the powers of the number of iterations. Thus this result can
3.2 CHEBYSHEV REGRESSION FOR DIAGONAL SHIFTS
For the shifted setting, we can use the same approach to prove that E,SSOR + cIn, b, w) is Lips- chitz w.r.t. the diagonal offset (c.f. Corollary F22); for n = 0(1) this implies regret Ö(T3/a VT) for the same discretization-batt based algorithm as Section 2.4. While optimal for Lipschitz functions, the method does not readily adapt to nice data, leading to various smoothed comparators (Krishnamurthy et al., 2019; Majzoubi et al., 2020; Zhu & Mineiro, 2022): however, as we wish to compete with the true optimal policy, we stay in the original setting and instead highlight how this section's semi-stochastic analysis allows us to study a very different class of bandit algorithms.
In particular, since we are now working directly with the cost function rather than an upper bound, we are able to utilize a more practical regression-oracle algorithm, SquareCB (Foster & Rakhlin, 2020). It assumes a class of regressors h |Cmin Cmin + C1 x [d] H [0, 1] with at least one function that perfectly predicts the expected performance S.SSOO + cI,, b. gig) of each action gii given the context c; a small amount of model misspecification is allowed. If there exists an online algorithm that can obtain low regret w.r.t. this function class, then SquareCB can obtain low regret w.r.t. any policy.
8

---

Algorithm 3: ChebCB: SquareCB with a follow-the- oracle and polynomial regressor class.
Input: solver SOLVE: R x R xΩ Z,00. instance sequence {(A,, b, b,)) IR"x R", context sequence {ct) 1=11 C [Cmin Cmin + C₁, learning rate yi > 0, parameter grid g E S24 Chebyshev polynomial features f: [(omin Cmin +C] mm+++, normalizations K, L,, N >
for f == ,,..., T' do
0: + arg min 18 10 KND
(0,
Re RN
Vield [d]
// update models
S f(c)) ViE [d]
11 compute model predictions
arg (e St1]
PM f 2+n($1y 8(,1)) Vi
1/ compute probability of each action
Pi+)] + 1 - L igti Pisl
// sample action sample is E [d] w.p. Pis and set wh -iae // run solver and update cost
hi t SOLVE A,, by, w))
To apply it we must specify a suitable class of regressors, bound its approximation error, and specify an algorithm attaining low regret over this class. Since terms of the Chebyshev series suffice to ap- proximate Lipschitz function with error (l/m)), we use Chebyshev polynomials in with learned coefficients-i.e. models (0, f(c) = En= Py)), where P3 is the jth Chebyshev polynomial- as our regressors for each action. To keep predictions bounded, we add constraints |راا| = O(1/5), which we can do without losing approximation power due to the decay of Chebyshev series coef- ficients. This allows us to show O(dm log regret for Hello-Thh--eader via Hazan et al. (2007, Theorem 5) and then apply Foster & Rakhlin (2020, Theorem 5) obtain the following guarantee: Theorem 3.2 (Corollary of Theorem C.4). Suppose Cmin -Amin (A). Then Algorithm 3 with appropriate parameters has regret w.E.t any policy f
[[mmin Cain + +C| + S2
dv Tvnn m
E
(W)))
Σ (f So dv ddvmnT++
Setting d = (T2/LL) and m= (TT/111 yields O(T9/11 regret, so asymptotically attain instance-optimal performance, albeit at a rather slow rate. The rate n is also worse than e.g. our semi-stochastic result for comparing to fixed w (c.f. Theorem 3.1), although to obtain this the latter algorithm uses d = O vn) grid points, making its overhead nontrivial. We compare ChebCB to the Section 2.4 algorithm based on Tsallis-INF (among other methods), and find that, despite the former's 's worse guarantees, it seems able to converge to an instance-opting policy much faster than the latter.
CONCLUSION AND LIMITATIONS
We have shown that bandit algorithms provably learn to parameterize SOR, an iterative linear system solver, and do as well asymptotically as the best fixed w terms of either (a) near-asymptotic measure of cost or (b) expected cost. We further show that modern contextual bandit method attains suariintaace-pptmaa performance. Both procedures require only the iteration count as feedback and have limited computational overhead settings, making them practical to deploy. Furthermore, the theoretical ideas in this work especially the use of contextual bandits for taking advantage of instance structure and Section 3. conversion of anti-concentrated Lipschitz criteria to Lipschitz expected costs-have the strong potential to be applicable to other domains of data-driven algorithm design.
At the same time, only the near-asymptotic results yield reasonable bound on the instances needed to attain good performance, with the rest having large spectral and dimension-dependen factors; the latter is the most obvious area for improvement. Furthermore, the near-asymptotic upper bounds are somewhat loose for sub-optimal w and for preconditioned CG, and as discussed in Section 2.4 do not seem amenable to regression- based CB. Beyond this, a natural direction is to attain semi-stochastic results for non-stationary solvers like preconditioned CG, or either type of result for the many other al- gorithms in scientific computing. Practically speaking, work on multiple parameters- e.g. the spectral bounds used for Chebyshev semi-iteration, multiple relaxation parameters for Block- Soock--OR-would likely be most useful. A final direction is to design online learning algorithms that exploit properties of the losses beyond I. Lipschitzness, or CB algorithms that take better advantage of such functions.

---

ACKNOWLEDGMENTS
We thank Akshay Krishnamurthy and Ainesh Bakshi for helpful feedback. This work was supported in part by National Science Foundation grants IIS-1705121, IIS-1838017, IIS-1901403, IIS-2046613, 12471, and OAC-2203821, the Defense Advanced Research Projects Agency under cooperative agreement HR001 12020003, TCS Presidential Fellowship, and funding from Meta, Morgan Stanley, Amazon, Google, and Jane Street. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of any of these funding agencies.
REFERENCES
Jacob Abernethy, Chansoo Lee, and Ambuj Tewari. Fighting bandits with a new kind of smoothness. In Advances in Neural Information Processing Systems, 2015.
Brandon Amos. Tutorial on amortized optimization. Foundations and Trends in Machine Learning, 16(5):592-732, 2023.
Hartwig Anzt, Edmond Chow, Jens Saak, and Jack Dongarra. Updating incomplete factorization preconditioners for model order reduction. Numerical Algorithms, 73:61 1-630, 2016.
Sohei Arisaka and Qianxiao Li. Principled acceleration of iterative numerical methods using machine learning. Proceedings of the 40th International Conference on Machine Learning, 2023.
Peter Auer, Nicolò Cesa-Bianchi, Yoav Freund, and Robert E. Schapire. The nonstochastic multiarmed bandit problem. SIAM Journal of Computing, 32:48-77, 2002.
Owe Axelsson. Iterative Solution Methods. Cambridge University Press, 1994.
Maria-Florina Balcan. Data-driven algorithm design. In Tim Roughgarden (ed.), Beyond the Worst- Case Analysis of Algorithmss Cambridge University Press, Cambridge, UK, 2021.
Maria-Florina Balcan, Travis Dick, and Ellen Vitercik. Dispersion for data-driven algorithm design, online learning, and private optimization. In 59th Annual Symposium on Foundations of Computer Science, 2018
Maria-Florina Balcan, Travis Dick, and Wesley Pegden. setting. Proceedings of the Conference on Uncertainty Artificial Intelligence, 2020.
optimization in the dispersed
Maria-Florina Balcan, Dan DeBlasio, Travis Dick, Carl Kingsford, Tuomas Sandholm, and Ellen Vitercik. How much data is sufficient to learn high-performing algorithms? Generalization guarantees for data-driven algorithm design. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computingg. 2021.
Maria-Florina Balcan, Mikhail Khodak, Dravyansh Sharma, and Ameet Talwalkar. Provably tuning the ElasticNet across instances. In Advances in Neural Information Processing Systems, 2022.
Peter Bartlett, Piotr Indyk, and Tal Wagner. Generalization bounds for data-driven numerical linear algebra. In Proceedings of the 35th Annual Conference on Learning Theory, 2022.
Manuel Baumann and Martin B. van Gijzen. Nested Krylov methods for shifted linear systems. SIAM Journal on Scientific Computing, 37:890-$111 2015.
Stefania Bellavia, Valentina De Simone, Daniela di Serafina, and Benedetta Morini. Efficient preconditioner updates for shifted linear systems. SIAM Journal on Scientific Computing, 33: 1785-1809, 2011.
Alina Beygelzimer, John Langford, Lihong Li, Lev Reyzin, and Robert E. Schapire. Contextual bandit algorithms with supervised learning guarantees. In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, 2011.
Sébastien Bubeck and Nicolò Cesa- Bosa-Bianchi. Regret analysis of stochastic and nonstochastic multi- armed bandit problems. Foundations and Trends in Machine Learning, 5(1):1-122, 2012.
10

---

Nicolò Desa-Bianchi and Gábor Lugosi. Prediction, Learning. and Games. Cambridge University Press, 2006.
Justin Y. Chen, Sandeep Silwal, Ali Vakilian, and Fred Zhang. Faster fundamental graph algorithms via learned predictions. In Proceedings of the 40th International Conference on Machine Learning, 2022.
Xinyi Chen and Elad Hazan. A nonstochastic control approach to optimization. arXiv, 2023.
Elliott Ward Cheney. Introduction to Approximation Theory. Chelsea Publishing Company, 1982.
Giulia Denevi, Carlo Ciliberto, Riccardo Grazzi, and Massimiliano Pontil. Learning-to- coaiing-to-leann learn stochas- tic gradient descent with biased regularization. In Proceedings of the 36th International Conference on Machine Learning, 2019.
Michael Dinitz, Sungjin Im, Thomas Lavastida, Benjamin Moseley, and Sergei Vassilvitskii. Faster matchings via learned duals. In Advances in Neural Information Processing Systems, 2021.
Paul Dütring, Guru Guruganesh, Jon Schneider, and Joshua R. Wang. Optimal no-regret learning for one-sided lipschitz functions. In Proceedings of the 40th International Conference Machine Learning, 2023.
Louis W. Ehrlich. An ad hoc SOR method. Journal of Computational Physics, 44:31-45, 1981.
Lakhdar Elbouyahyaoui, Mohammed Heyouni, Azita Tajaddini, and Farid Saberi-Movaher On restarted and deflated block FOM and GMRES methods for sequences of shifted linear systems. Numerical Algorithms, 87:1257-1299, 2021.
Dylan J. Foster and Alexander Rakhlin. Beyond UCB: Optimal and efficient contextual bandits with regression oracles. In Proceedings of the 37th International Conference on Machine Learning, 2020.
Isaac Fried and Jim Metzler. SOR vs. conjugate gradients in finite element discretization. Interna- tional Journal for Numerical Methods in Engineering, 12:1329-133 1332, 1978.
Andreas Frommer and Uew Glässner. Restarted GMRES for shifted linear systems. SIAM Journal on Scientific Computing, 19:15-26, 1998.
Gene H. Golub and Qiang Ye. Inexact preconditioned conjugate gradient method with inner-outer iteration. SIAM Journal on Scientific Computing, 21:1305-133 1320, 1999.
Anne Greenbaum. Iterative Methods for Solving Linear Systems. Society for Industrial and Applied Mathematics, 1997.
Rishi Gupta and Timothy Roughgarden. A PAC approach to application-specific algorithm selection. SIAM Journal on Computing, 46(3):992-1017, 2017.
Wolfgang Hackbusch. Iterative Solution of Large Sparse Systems of Equations. Springer International Publishing, 2016.
Elad Hazan, Amit Agarwal, and Satyen Kale. Logarithmic regret algorithms for online convex optimization. Machine Learning, 69:169-192, 2007.
George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 11422-440 2021.
Michael Kazhdan, Jake Solomon, and Mirela Ben-Chen. Can mean-curvature flow be modified to be non-singular? In Eurographics Symposium on Geometry Processing 2012, 2012.
Mikhail Khodak, Maria-Florina Balcan, and Ameet Talwalkar. Adaptive gradient- meta-learning methods. In Advances in Neural Information Processing Systems, 2019.
Mikhail Khodak,
Balcan, Ameet Talwalkar, and Sergei Vassilvitskii. Learning pre- dictions for algorithms with predictions. In Advances in Neural Information Processing Systems, 2022.
11

---

John B. King, Samim Anghaie, and Henry M. Domanus. Comparative performance of the conjugate gradient and SOR methods for computational thermal hydraulics. In Proceedings of the Joint Meeting of the American Nuclear Society and the Atomic Industrial Forum, 1987.
Robert Kleinberg. Nearly tight bounds for the continuum-armed bandit problem. In Advances in Neural Information Processing Systems, 2004.
Akshay Krishnamurthy, John Langord, Alexandrs Slivkins, and Chicheng Zhang. Contextual ban- dits with continuous actions: Smoothing, zooming, and adapting. In Proceedings of the 32nd Conference on Learning Theory, 2019.
John Lafferty, Han Liu, and Larry Wasserman. https//wwww.ttt.cmu.dddd Barry/esmm/Conceentaatinnn 2010.
Statistical
machine learning.
Sören Laue, Matthias Mitterreiter, and Joachim Giesen. Computing higher order derivatives of matrix and tensor expressions. In Advances Neural Information Processing Systems, 2018.
Randall J. Ve Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM, 2007.
Yichen Li, Peter Yichen Chen, Tao Du, and Wojciech Matusik. Learning preconditioners for conjugate gradient PDE solvers. In Proceedings of the 40th International Conference on Machine Learning, 2023.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadeneshelii Burigede Liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In Proceedings of the 9th International Conference on Learning Representa- fions, 2021.
Tyler Lu, Dávid Pál, and Martin Pál. Contextual multi-armed bandits. In Proceedings of the 13th International Conference on Artificial Intelligence and Statistics, 2010.
Ilay Luz, Meirav Galun, Haggai Maron, Ronen Basri, and Irad Yavnch. Learning algebraic multigrid using graph neural networks. In Proceedings of the 37th International Conference Machine Learning, 2020.
Maryam Majzoubi, Chicheng Zhang, Rajan Chari, Akshay Krishnamurthy, John Langford, and Alexandrs Slivkins. Efficient contextual bandits with continuous actions. In Advances in Neural Information Processing Systems, 2020.
Donald Marquardt. An algorithm for least-s estimation of nonlinear parameters. SIAM Journal on Applied Mathematics, 1 11(2):31-441, 1963.
Tanya Marwah, Zachary C. Lipton, and Andrej Risteski. Parametric complexity bounds for approx- imating PDEs with neural networks. Advances in Neural Information Processing Systems, 2021.
Michael Mitzenmacher and Sergei Vassilvitskii. Algorithms with predictions. In Tim Roughgarden (ed.), Beyond the Worst-Case Analysis of Algorithms. Cambridge University Press, Cambridge, UK, 2021.
Cameron Musco and Christopher Musco. Randomized block Krylov methods for stronger and faster approximate singular value decomposition. In Advances in Neural Information Processing Systems, 2015.
Michael L. Parks, Eric de Sturler, Greg Mackey, Duane D. Johnson, and Spandan Maiti. Recycling Krylov subspaces for sequences of linear systems. SIAM Journal on Scientific Computing, 28(5): 1651-1674, 2006.
Shinsaku Sakaue and Taihei Oki. Diisreee-oomee--aallis-----e framework for warm-starting algorithms with predictions. In Advances in Neural Information Processing Systems, 2022.
Rajiv Sambharya, Georgina Hall, Brandon Amos, and Bartolomeo Stellato. End-to-end learning to warm-start for real-time quadratic optimization. In Proceedings of the 5th Annual Conference on Learning for Dynamics and Control, 2023.
12

---

Nikunj Saunshi, Yi Zhang, Mikhail Khodak, and Sanjeev Arora. A sample complexity separation be- tween non-convex and convex meta-learning. Proceedings of the 37th International Conference on Machine Learning. 2020.
David Simchi-Levi and Yunzong Xu. Bypassing the monster: faster and simpler optimal algorithm for contextual bandits under realizability. Mathemalics of Operationn Research, 47, 2021.
Ali Taghibakhshi, Scott MacLachlan, Luke Olson, and Matthew West. Optimization bpgimization-basedd algebraic multigrid coarsening using reinforcement learning. Advances in Neural Information Processing Systems, 2021.
Jurjen D. Tebbens and Miroslav Tûma. Efficient preconditioning of sequences of nonsymmetric linear systems. SIAM Journal on Scientific Computing. 29:1918-1941, 2007.
James William Thomas. Numerical Partial Differential Equations. Springer Science+Business Media, 1999.
Lloyd Trefethen. Gauss quadrature better than Clenshaw-Curtis? SIAM Review, 50(1):67-87, 2008.
Lloyd N. Trefethen and Mark Embree. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, 2005.
L. Dale Van Vleck and D. J. Dwyer. Successive overrelaxation, block iteration, and method of conjugate gradients for solving equations for multiple trait evaluation of sires. Jorunal of Dairy Science, 68:760-767, 1985.
Rui-Rui Wang. Qiang Niu, Xiao-Bin Tang, and Xiang Wang. Solving shifted linear systems with restarted GMRES augmented with error approximations. Computers & Mathematics with Applica- tions, 78: 1810-19119,8 2019.
Zbigniew I. Wožnicki. On numerical analysis of conjugate gradient method. Japan Journal of Industrial and Applied Mathematics, 10:487-519, 1993.
Zbigniew I. Wožnicki. On performance of SOR method for solving nonsymmetric linear systems. Journal of Computational and Applied Mathematics, 137:145-176, 2001.
David M. Young, Iterative Solution of Large Linear Systems. Academic Press, 1971.
Yinglun Zhu and Paul Mineiro. Contextual bandits with smooth regret: Efficient learning in contin- uous action spaces. In Proceedings of the 39th International Conference Machine Learning. 2022.
Julian Zimmert and Yevgeny Seldin. Tsallis-NNF An optimal algorithm for stochastic and adversarial bandits. Journal of Machine Learning Research, 22:1-49, 2021.
13

---

A RELATED WORK AND COMPARISONS
Our analysis falls mainly into the framework of data-driven algorithm design, which has a long history (Gupta & Roughgarden, 2017: Balcan, 2021). Closely related is the study by Gupta & Roughgarden (2017) of the sample complexity of learning the step-size of gradient descent, which can also be used to solve linear systems. While their sample complexity guarantee is logarithmic in the precision 1 1/s, directly applying their Lipschitz-like analysis in a bandit setting yields regret with a polynomial dependence; note that a typical setting of E is 10-a Mathematically, their analysis relies crucially on the iteration reducing error at every step, which is well-known not to be the case for SOR (e.g. Trefethen& & Embree (2005, Figure 25.6)). Data-driven numerical linear algebra was studied most explicitly by Bartlett et al. (2022), who provided sample complexity framework applicable to many algorithms; their focus is on the offline setting where an algorithm is learned from a batch of samples. While they do not consider linear systems directly, in Appendix A. we do compare the guarantee their framework implies for SOR; we obtain similar sample complexity with an efficient learning pro- cedure, at the cost of strong distributional assumption on the target vector. Note that generalization guarantees have been shown for convex quadratic bogrramming-whinn which we do not consider because for
subsumes linear systems-
by Sambharya al. (2023); they focus on learning-to-initee high precisions the initialization quality usually does not have strong impact on cost. Note that all of the above work also does not provide end-to-end guarantees, only e.g. sample complexity bounds.
Online learning guarantees were shown for the related problem of tuning regularized regression by Bal- can et al. (2022), albeit in the easier full information setting and with the target of reducing error rather than computation. Their approach relies on the dispersion technique (Balcan et al., 201 8). which often involves showing that discontinuitics in the cost are defined by bounded-degrce polynomials (Balcan et al., 2020). While possibly applicable in our setting, we suspect using it would lead to unacceptably high dependence on the dimension and precision, as the power of the polynomials defining our deci- sion boundaries is On lo Lastly, we believe our work is notable within this field as a first exam- ple of using contextual bandits, and in doing so competing with the provably instance-optimal policy,
Iterative (discrete) optimization has been studied in the related area of Saaming-augmenter algorithms (a.k.a. algorithms with predictions) (Dinitz et al., 2021; Chen et al., 2022; Sakaue & Oki, 2022), which shows data-dependent performance guarantees as a function of (learned) predictions (Mitzenmacher & Vassilvitskii, 2021); these can then be used as surrogate losses for learning (Khodak et al., 2022). Our construction of an upper bound under asymptotic convergence is inspired by this, although unlike previous work we do not assume access to the bound directly because it depends on hard-to-come to-comptte spectral properties. Algorithms with predictions often involve initializing a computation with a prediction of its outcome, e.g. vector near the solution A --;b; we do not consider this because the runtime of SOR and other solvers depends fairly weakly on the distance to the initialization.
A last theoretical area is that of gradient-based meta-learning, which studies how to initialize and tune other parameters of gradient descent and related methods (Khodak et al., 2019; Denevi et al., 2019; Saunshi et al., 2020; Chen & Hazan, 2023). This field focuses on learning- thhooretc notions of cost such as regret or statistical risk. Furthermore, their guarantees are usually on the error after fixed number of gradient steps rather than the number of iterations required to converge; targeting the former can be highly suboptimal in scientific computing applications (Arisaka & Li, 2023). This latter work, which connects meta-learning and data-driven scientific computing, analyzes specific case studies for accelerating numerical solvers, whereas we focus on a general learning guarantee.
Empirically, there are many learned solvers (Luz et al., 2020; Taghibakhshi et al., 2021; Li et al., 2023) and even full simulation replacements (Karniadakis et al., 2021; Li et al., 2021); to our knowledge, theoretical studies of the latter have focused on expressivity (Marwah al., 2021). Amortizing the cost on future simulations (Amos, 2023), these approaches use offline computation to train models that integrate directly with solvers or avoid solving linear systems altogether. In contrast, the methods we propose are online and lightweight, both computationally and in terms of implementation; unlike many deep learning approaches, the additional computation scales slowly with dimension and needs only black-box access to existing solvers. As result, our methods can be viewed as reasonable baselines, and we discuss an indirect comparison with the (G-pprcondittione--lammi et al. (2023) in Appendix G. Finally, note that improving the performance of linear solvers across sequence of related instances has seen a lot of study in the scientific computing literature (Parks et al., 2006; Tebbens & Túma, 2007; Elbouyahyaoui et al., 2021). To our knowledge, this work does not give explicit guarantees on the number of iterations, and so a direct theoretical comparison is challenging.
approach of Li
14

---

A.1 SAMPLE COMPLEXITY AND COMPARISON WITH THE GOLDBERG-JERRUNE FRAMEWORK
While not the focus of our work, we briefly note the generalization implications of our semi-stochastic 3) i.i.d. samples from a distribution D over matrices A, satisfying the assumptions in Section 2.1 and truncated Gaussian targets b arg minges Σ SSOR(A, be, over uniform grid [1, [ of size d = Ō(VnT) will be a-suboptimal w.p. 1- 5:
analysis. Suppose for any a > O we have T = O oz polylog
Then empirical risk minimization
Corollary A.1. Let D be a distribution over matrix-vector pairs (A.b) E Ruxa x RR where A satisfies the SOR conditions and for every A the conditional distribution of D given A over B" is the truncated Gaussian. For every T > consider the algorithm that draws T samples
(A,, bi) D and outputs @ arg mingeg Et= (S)B;(w), where g(i] = =1+ (Wmax 1) 2
and d LVT for I as in Corollary F.I. Then T = Ô po2poly 28 sumples suffice to ensure
K
KK
ED SSOR A, b, S minwe(l [[max] SSOR A, b, + a holds w.p. 2 1 8.
Proof. A standard covering bound (see e.g. Lafferty et al. (2010, Theorem 7.82)) followed by an application of Corollary F.1 implies that w.p. 2
2d
EpSsOR(A, z min EpSSOR(SSOAAA A, b, + 3K log
Weg
T
8
2
2d
min BA (,, +K V T WEg
log
8
2
2d
min
EA EbSSORRA, b, +
+ 3K
T
log
(10)
I. d
2
2d
min
EpSSOR(A, b,
3K
log
T
2
2LT
< min [1(1.
EpSSOR(A, b, b, +4K
log
T
K8
Noting that by Corollary F.1 we have L = O Vn) = OlV/nol' ?) yields the result.
L
This matches directly applying the GJ framework of Bartlett al. (2022, Theorem 3.3) to our problem:
Corollary A.2. In the same setting as Corollary A.I but generalizing the distribution to any one whose target vector support is vn-bouuded. empirical risk minimization (running a =
arg nigmin (1mmmma)) Σt (SOR((w) has sample complexity Ô poooyoo
Proof. For every (A, b) pair in the support of D and any R it is straightforward define GJ algorithm (Bartlett et al., 2022, Definition 3.1) that checks if SSOR(A, b, > r by computing [|rx(w)]13 = C+ b 1/3_ -/ degree 2k folynomial-for every k: < [r] and returning "True" if one of them satisfics | (w)]l < 2 and "False" otherwise (and automatically return "True" for r 2 K and "False" for r < 1). Since the degree of this algorithm is at most 2K, the predicate complexity is at most K, and the parameter size is 1, by Bartlett et al. (2022, Theorem 3.3) the pseudodimension of {SSOR(: w) w E [1. (mxx)) }} is (log K). Using the bounded assumption on the target vector -completes the proof.
SSCR < K = (log 1
At the same, recent generalization guarantees for tuning regularization parameters of linear regression by Balcan et al. (2022, Theorem w.o applied dual function analysis (Balcan et al., 202e>-have a quadratic dependence on the instance dimension. Unlike both results-which use uniform convergence-out bound also uses (theoretically) efficient learning procedure, at the cost of a strong (but in our view reasonable) distributional assumption on the target vectors.
15

---

A.2 APPROXIMATING THE SPECTRAL RADIUS THE JACOBI ITERATION MATRIX
Because the asymptotically optimal w is a function of the spectral radius 3 - p(Mi) ) of the Jacobi iteration matrix, reasonable baseline is to simply approximate 3 using eigenvalue solver and then run SOR with the corresponding approximately best w. is difficult to compare our results to this approach directly, since the baseline will always run extra matrix iterations while bandit algorithms will asymptotically run no more than the comparator. Furthermore, M1 is not a normal matrix, class for which turns out to be surprisingly difficult find bounds on the number of iterations required to approximate its largest eigenvalue within some tolerance a > 0.
comparison can be made in the diagonal offset setting by modifying this baseline somewhat and making the assumption that A has constant diagonal, so that M1 is symmetric and we can use randomized block- Krylov to obtain 3 satisfying 1 - 3 Oss) in ((1/ VE) iterations w.h.p. (Musco & Musco, 2015, Theorem 1). To modify the baseline, we consider preprocessing algorithm which discretizes [Cmin= Cmin + C] into d grid points, runs & iterations of randomized block- Krylov on the Jacobi iteration matrix of each matrix + cI, corresponding to offsets c in this grid, and then for each new offset C₃ set wt using the optimal parameter implied by the approximate spectral radius of the Jacobi iteration matrix of A + cI,, corresponding to the closest e in the grid. This algorithm thus does (dke) matrix-vecton products of preprocessing. and since the upper bounds U, are H-Hölder w.r.t. while the optimal policy is Lipschitz w.r.t. 8 over an appropriate domain [1, wiax]] will w.h.p. use at most Ó V1/12 + 1/d) more iterations at each step /. E [T] compared the optimal policy. Thus w.h.p. the total regret compared to the optimal policy w"
T Σ
T
(w)) 7 (dk (dk++ T/l++T/+T// + (vllw(a))
(11)
Setting d VT and k VT yields the rate Ö(T3/4), which can be compared directly to O(T3/4) rate for the discretized Tsallis-INF algorithm in Theorem 2.2. The rate of approximating p(M)) thus matches that of our simplest approach, although unlike the latter (and also unlike ChebCB) it does not guarantee performance as good the optimal policy in the semi-stochastic setting, where w' might not be optimal. Intuitively, the randomized block- block-Krylov baseline will also suffer from spending computation on points [Cmins Cmin C* that does not end up seeing.
16

---

Algorithm 4: General form of Tsallis-INF, The probabilities can be computed using Newton's method (e.g. Zimmert& & Seldin (2021, Algorithm 2)).
Input: loss sequence (le : (a, bj +> [0, K K]) action set g E (a, bjd, step-sizes ",,..., Or // initialize vector of cumulative losses
k <- 0a
for f ,,..., do
< arg m (k. - 1K Ex=
T/ ii=1 PP
// compute probabilities // sample index of an action // play action and update losses
sample [d] with probability P(t]
kid kppl (((())))P(.))
B SEMI-LIPSCHITZ BANDITS
We consider a sequence of adaptively chosen loss functions l1...., [a, + [0, K] on an interval [a. b] C IR and upper bounds u1...., ur : Ja, b] + R satisfying u, (x)x) 2 (x) A E [T. [a, bj, where [T] denotes the set of integers from 1 to T. Our analysis will focus on the Tsallis-INF algorithm of Abernethy al. (2015), which we write in its general form in Algorithm 4, although the analysis extends easily to the better-known (but sub-optimal) Exp3 (Auer et al., 2002). For Tsallis-NNF, the following two facts follow directly from known results:
Theorem B.1 (Corollary of Abernethy al. (2015, Corollary 3.2)). 1f 1/VTVt [T] then
Algorithm 4 has regret E Σ (((g)))) min;g[4] Σ & (g:)) S 2Ky2dT.
Theorem B.2 (Corollary Zimmert& Seldin (2021, Theorem I)). If Th = 2/ E [T] then
Algorithm 4 has regret E Σí= ((()))) min:e[4] 2 (((a) < 4K /dT +
We now define generalization of the Lipschitzness condition that trivially generalizes regular L-LLipschizz functions, as well as the notion of one-sided Lipschitz functions studied in the stochastic setting by Dütting et al. (2023).
Definition B.1. Given a constant 2 and a point & E [a, b). we say a function f f: la, + R is (L,2) se..)))smm-Liiphhizz f(x) f(v) S D|x yi v x, y s.t. Jx 2
2/-
We now show that Tsallis INF with bandit access to l, on discretization of [a, b] attains ((72/3)
regret w.r.t. any fixed æ E [a, bj evaluated by any comparator sequence of semi-Lipschitz upper
bounds uy. Note that guarantees for the standard comparator can be recovered by just setting 4 = ue V E [T]. and that the rate is optimal by Kleinberg (2004, Theorem 4.2).
Theorem B.3. u 2 l is (L,,2) semis-hiih V t [T] then Algorithm 4 using action space E
[a, bjd s.t. g[ = a + biiiii dd 1] and gid = z has regret
T
b
T
E
Σ ((g())
min xE(a,)]
ux(x) S 2K 2dT
2dT
ΣL.
Lr
(12)
d
Setting d = 3/ (b (b-a)P1 P.PT for T. 1 Et T4 yields the bound 3 2/2(1 TTTKKT22
2K
Proof. Let [- denote rounding to the closest element of g in the direction of z. Then for x e [a,b]
we have |[ [ 2| $ |x
2/ and |[>|g x
S b-a d so applying Theorem B.1 and this fact yields
T
E
((gRR)) 2KK2dT min 1
& (g I 2K
2dT
min iE|d
ue (g(4)
2K 22 + min zE(a,b)
u [x: g)
b
T
T
S 2K 22KK2dl +
Σ
L, +
min
u4 (x)
d
[a
(13)
17

---

Algorithm 5: Contextual bandit algorithm using instances of Tsallis-NN over grid of contexts.
Input: loss sequence (le (a, + [0, K K] context sequence (4 C [c,c C). action set g (a. b)s, discretization h E [,,++ CI
for j = ... m do
(gF ( // start m instances of Algorithm Aj Tsallis INF
for / 1....; do
Jt = min arg min je(ml |b al // pick element of h closest to 9
+ Aje
// get action from jthh instance of Algorithm // pass loss to jihh instance of Algorithm
((K.)) + Aj
For contextual bandits, we restrict to (I,, semissmimmipiiiizzz functions and T.p-Lipsc policies, obtaining ((T3/4) regret; this rate matches known upper and lower bounds for the case where losses Lipschitz in both actions and contexts (Lu et al., 2010, Theorem 1), although this does not imply optimality of our result.
Theorem B.4. If Ut 2 Ee is (Lr, siemi--iiii and & E lc, CVIC [T] then Algorithm using action space gis] = + J and h11 = c+ ci 1) as the grid of contexts has regret w.r.t.
any Lipscchtt policy f [c c CJ b0,
T
[CLy
b
T
E
(g(a))
Σ ux(T(e))) & m + 4K dmmT +
(14)
d
1=1
(b-
CIITT 4(baajk2 yields regret 4 4/4K a ACb,K2 K2 m V
1
Setting d
4(b-a)k2.
Proof. Define I be the operation of rounding to the closest element of h, breaking ties arbitrarily. and set [T], = {t E [T] : [cel = hp)]. Furthermore, define [x] to be the smallest element g[i] in CL2
g s.t. x 2m
S[i] (or max4e[d] Big if such an element does not exist).
T
E
((gg))
E
Σ
(gg..)
min sEld
((g
+ min Σ (gg())
ieid
j=1 tE|T
LE[TT
1E(T]
m
K / / + min
Σ (gg)
m
iEid]
1=1
rE[T]
(15)
m + 4K ddmT
min EE
u((g))
f=1
(E(T];
m 4K <V dmmT
u ([f((c[[)))))
where the first inequality follows by Theorem B.2, the second applies Jensen's inequality to the left term and u, 2 l, on the right, and the last uses optimality of each for each j. Now since ff is Lipischhi we have by definition of [- that (f(c)) - f([caJn) < CL2 This in turn implies that f(a)) < [f([c]])] < f(c:) + CLI b-a by definition of g and Flg. Since u, is (L,, b) siii---iiiss the result follows.
2m
18

---

C CHEBYSHEV REGRESSION FOR CONTEXTUAL BANDITS
C.I PRELIMINARIES
We first state Lipschitz approximation result that is standard but difficult-tofindd formally. For all E20 we will use (x) = cos(j arccos(z)) to denote the jth Chebyshev polynomial of the first kind. Theorem C.1. Let f |=1] + [ be a K K-bounded, L-Lipschitz function.. Then for each integer m 2 0 there exists 0 E R=+ satisfying the following properties:
1. [00] S K and [0/]] - 2L/jV3 j E [m]]
2. maxze]+11 ((x)) - -Ø - D1 P3( I
タ+ 1og(2m++1) L
m+1
Proof. Define Quj = S1 f(z) de and for each j E [m]] let 0b = - f())((()) da be the jth -I-
vI-z
Chebyshev coefficient. Since 1 VI-x Theorem 4.2) we also have
de
trivially have |0|o]| K and by Trefethen (2008,
2
[f()))
2L dr 2L/j Tj
dx -1-
(16)
Tj.
11-12
for all j [m]. This shows the first property. For the second, by Trefethen (2008, Theorem 4.4) we have that
m Σ
4log(2m + +1)
max xE|-1,11
f(x)
Py())
$ 2+
max xE[=1
[f(x) Pm (x))
x²
j=0
log (2m + 1)
La 2 + 1)
log (2m + 1)
$ 2
L
π
+
(17)
where Pm is the (at most) m-degree algebraic polynomial that best approximates on [+1] and the second inequality is Jackson's theorem (Cheney, 1982, page 147).
Corollary C.1. Let f [a, [+]] be a K-bounded, L-Lipschitz function the interval [a, b]. Then for each integer m 2 0 there exists 0 E Rm+11 satisfying the following properties:
b]
1. 00/] S K and 10 L(b-a)
2. maxzela (xx) - -5 O1n Pr (1-22220000 11/
s+ 3 (2)
L(b a)
2(m+l)
(1 1)+ a). that g [+1] [+K] is K-bounded and L Lez. Lipschitz. Applying Theorem C.1 yields the result.
Proof. Define g(x)
We next state regret guarantees for the SquareCB algorithm of Foster & Rakhlin (2020) in the non-realizable setting:
Theorem C.2 (Foster & Rakhlin (2020, Theorem 5)). Suppose for any sequence of actions 01, 01...., or an online regression oracle A playing regressors hi,,..,,lv lur H has regrel guarantee
T
RT 2 (lelee)))) ax)))
min KCH
(El(c, az) - h(c,, az))
(18)
If all losses and regressors have range [0. and 3 h E H x.l. Bl((a) Mcc, a) + onlee, a) for [ay(a)) S then Algorithm 6 with learning rate 2 dl / 2a²T" has expected regret w.r.t the the optimal policy f |a, bj + g bounded as
E
e (g(.)) - Σ ((((())) s 2 5d
(19)
19

---

Algorithm 6: SquareCB method for contextual bandits using online regression oracle.
Input: loss sequence (le > [0, 10] context sequence (ca [c.)21, learning rate >> O, online regression oracle A for f 1...., T do
Sig f A(c,, g(a) ViE [d]
// compute oracle prediction
i* f arg min, Stg
PK -- .p1) Vi
// compute action probabilities
Pic] 4-
jp' P
sample [d] with probability P(a] // sample index of grid point ((c. - A // pass context, action, and loss to oracle
Algorithm SquareCB method for Lipschitz contextual bandits using Follow-the-Leader
Input: loss sequence {6e : [a, b] + [0, K]]== context sequence {er E [c, + CT) learning rate > 0, action set g E [a, bl", featurizer f. : cc]c+] + IR",, normalizations L, N > 0 for ,,..,, do
Σ (10,((c)))
(( 2 KN
0. + arg min
E // update models
p0 1o RNy sE(-11)
S (0,, f(c)) Vie [d]
// compute model predictions
arg min, SK)]
P f d+ y )p.))
Vi
// compute action probabilities
Ple] ii PM
sample i E [d] with probability P(a] and play action g[ú]]
SquareCB requires an online regression oracle to implement, for which we will use the Follow-the- Leader scheme. It has the following guarantee for squared losses:
Theorem C.3 (Corollary of Hazan et al. (2007, Theorem 5)). Consider the follow-the-leader al-
gorithm, which sequentially sees feature-target pairs (x1, v1 y1), ... (XT, yT) E x [0, 1] for some subset X C. [0, 1l" and at each step sets 0,+1 = arg migmingee Ei=(( 0) - yy)) for some subset R". This algorithm has regret
T
XDT
Σ ((x,, 0,) - yz)² - -in2-min((6 (xx, 0) yx)) < 4B2n 1 log
(20)
d58
2B
for De the diameter maxg,8 110 12 O, X max.e(1] ||x||22 and max (xx, 0)}-
C.2 REGRET OF CHEBCB
Theorem C.4. Suppose El((x) is an LipLicschitz function actions [a, and an L.- Lipschitz function of contexts e E [c.c + C}. Then Algorithm 7 run with learning rate = 2/ /RT + 2a2T) for Rr and as in Equations 22 and 23, respectively, action set g( (b a) 1112 Chebyshev features (c)) = (ct,,,, and normalizations L L
and = -2+ ACL, (1 + log has regret w.r.t. any policy f lc,c C1++ [a, of K
L.TVd T.gT Σ (g.))) ((f(c))) = Ô [..dVmm
E
+
(21)
m
d
Setting d = 1(T2/11 and m= = (TT311)) yields regret Ó(max(L T T99111)
20

---

Proof. Observe that the above algorithm is equivalent to running Algorithm 6 with the follow-the-
leader oracle over d(m + -)-dimensional space & with diameter V N 1 4C
ak /K3+2UC-1222/3 features bounded by V
1 1+ Σj= Pf(c) & vm + and predictions bounded
KN
by (((()),)] S |10 ||))))) S N + most
2CLe KN
Σド= S 2. Thus by Theorem C.3 the oracle has regret
Rr = = 1) 1
T' d( 1)(K², 2C"
(22)
KN
Note that, to ensure the regressors and losses have range in [0, 1] we can define the former as (e, 8i4) = (f(c), 0,) + and the latter as KN 2 and Algorithm 7 remains the same. Furthermore, the error of the regression approximation is then
f +
+ 2 log(2m 1)
CL.
(23)
2KN (m + 1)
We conclude by applying Theorem C.2, unnormalizing by multiplying the resulting regret by KN, and adding the approximation error (b( due to the discretization of the action space. L
2d
21

---

D SOR PRELIMINARIES
We will use the following notation:
M I ID/=+L)-1 is the matrix of the first normal form (Hackbusch, 2016, 2.2.1) the matrix of the third normal form (Hackbusch, 2016, Section 2.2.3) C L, AA(DDu L) L) I, AW AM,.A is the defect reduction matrix
w. D/w
M. = In 2 DD/u +L=) +D + - AA the matrix of the first normal form for SSOR (Hackbuseh, 2016, 2.2.1)
W., -ー (D/W~ L)D- + L²² ) is the matrix of the third normal form for SSOR (Hackbusch, 2016, Section 2.2.3)
Cw = I, 2- A(D/w + (L")-"D(D/aa + + - I AW. AM is the defect reduction matrix for SSOR
|| ||2 denotes the Fuclidean norm of a vector and the spectral norm of a matrix
<(A) |A| A
A
ll2 denotes the condition number of a matrix A
p(X) denotes the spectral radius of matrix X
| = || ||A+x|| x A >
denotes the energy norm of a vector I" associated with the matrix
| | - ||2 denotes the energy of a matrix X Raxx associated with the matrix A > 0
We further derive bounds on the number of iterations for SOR and SSOR using the following energy norm estimate:
Theorem D.1 (Corollary of Hackbusch (2016, Theorem 3.44 & Corollary 3.45)). If A > 0 and where = =1 (L+
2
E (0, 2) then [| <
Corollary D.1. Let Kw be the maximum number of iterations that SOR needs to reach error > loog 2Ve(A)
Then for any E (0,2) we have K I 1+ upper bound in Theorem D.I.
where (A(A) is the square root of the log
Proof. By Hackbusch (2016, Equations 2.22c & B.28b) we have at each iteration is of Algorithm that
| ll
2 |A
| r S 22/x(A))
(24)
|roll
||rol12
Setting the r.h.s. equal e and solving for k yields the result.
Corollary D.2. Let Ka the maximum number of iterations that SSOR (Algorithm 8) needs to reach (absolute) error > O. Then for any E (0,2) we huve K2, < 1 + Vol (A) is the square root of the upper bound Theorem D.I.
log 21bl where log (A
Proof. By Hackbusch (2016, Equations 2.22c & B.28b) have each iteration of Algorithm that
|| ll S 2|| | |A A
rolle 2|]blla
(25)
Setting the r.h.s. equal to and solving for & yields the result.
22

---

0.8
0.8
V(1
1/e2
0.2
-0.1
0.1
0.2
0.3
0.4
Figure 3: Values of r and 33 for A + cI, for different c.
E
NEAR-ASYMPTOTIC PROOFS
E.1 PROOF OF LEMMA 2.1 AND ASSOCIATED DISCUSSION
Proof. For the first claim, suppose
min |
x> U(w). Then
bilaa blla
1ll2
E
s (p(C.) + + <E
(26)
|bll
|
so by contradiction we must have
min bcbbla<<llbll
S U(w). Now note that by Hackbuseh (2016,
Theorem 4.27) and similarity of C., and M., we have that p(C..) (wB + ww232 4(w
for 1
1-V and 1 otherwise. Therefore
wn have p Cw) s
p(C1) = 32 and w w have p(Cw) Whax 1. This concludes the second part of the first claim. The first part of the second claim follows because p (Cw) is decreasing on w w" For Since
the second part, compute the derivative [a,((w)) (-+(1-)(w-1) (o (1+(1111)) --- (w--))
(r-1) log
T (w 1) 2 by assumption-eiteee either by nonnegativity of w if r 2
or because
otherwise B² 2 (1 - implies + (1 (w 1) 2 (1 1/)
the
I+ 1+V1-33
derivative is increasing in w and so is most -[1-T) log" a
Note that the fourth item's restriction on 7 and 3 does not really restrict the matrices our analysis applicable to,
always re-define in Assumption 2.1 be least 1/e2. Our analysis does not strongly depend on this restriction; it largely done for simplicity and because it does not exclude too many settings interest, In-aartic,lar, we find for 10-8 that typically indeed larger than 1., and furthermore is likely quite high whenever / is small, as it suggests the matrix is near-diagonal and so w near one will converge extremely quickly (c.f. Figure 3).
E.2 PROOF OF THEOREM 2.1
Proof. The first bound follows from Theorem B.3 by noting that Lemma 2.1 implies that the functions
Wmas -semi-Lipschitz over [1, wmxx] and the functions SOR, - 1 < U, - lo -bounedd. To extend the comparator domain (0, (mac), note that Lemma 2.1.2 implies that all are decreasing on w E (0, 1). To extend the comparator domain again in the second bound, note that the setting of Waax implies that the minimizer 1 + 88 (1 + /1- -87) of each U, is at most (Wimxx 2).
U, - are (-(1-)) loge
log
logos or,
Wmory and all functions U, increasing
23

---

E.3 APPROXIMATING THE OPTIMAL POLICY
Lemma E.I.
Define A - A + for all E [Cinns 20), where Cmin Amii (A). Then w" (c)
1+
Be
is
60max (1 (/mxa)/ 1--0ixxx
(Avia (D) +min
-Lipschitz, where Beax max, Be
(I-
(Vi- +1)
the maximum over 3e = p(In -D (A+c)))(AAA
Proof. We first compute
0.3, = dcp(In - + (
Dedmax (I, (D+
-viO,([,, (D +))1]
(D++
2 (00 2c(D + cl,)
I. ( (A+ cI,)
- ?I, (D +c,))))-
2 2V((D+ (4 ivi
2 d c., (I,) iv ov cI)) 'v1
Vi (I, +D++ (D+ + c + cI,) - - + (D + cl,)-")v
cI.)
2
Df 3cL cI, iv
vf -3 cI.) cI. vv
(27)
The first component is positive and the matrix has eigenvalues bounded by (1+ > (D)+c
1+
while the last term is negative and the matrix has c 2 0 the positive component has spectral radius the middle term is negative, subtracting 2A from middle matrix shows that its magnitude is bounded by [1 + Be). the middle term positive- wwhic can
most 1+
-3e/2 only happen for negative c-its magnitude is bounded by (D)
3Amin (Di/2 Combining all Amia (D)+c*
terms yields bound 3 Amin (D)+c+1
(1 (1+ Bc). We then have that
(D)+
10.w" (c)l = 23./ / 32
2 10 = 63. 3e) 3.) 32
(Amia MD)
(28)
A (D) + (D)
(VI - 32
(V1- (
The result follows because 2x(1+x)/ v-1" decreases monotonically in
increases monotonically on r E [0, and the bound itself
(V1 (VI-22+1)*
24

---

S-shaped domai:: size 100, otfasie0 port tound
R-shaped domain: siep-100, offeet48
B-shaped domain: sizu-103, ottasteg
R-ahaped domain: siee-200, offeet-0.
fooed
tound)
foune
L-shapedd domales slea-75, offscte)
offeet40.5
L-shaped sloeo675
Lyshappd slesed75, offsese8.
con
Figure 4: Comparison of actual cost of running SSOR-preeondittioned CG and the upper bounds computed in Section E.4 as functions of the tuning parameter w E [2v2 - 1.9] on various domains.
Theorem E.1. Suppose Ot E [Cmins Cmin + V [T], where Cmin > -Amin (A), and define Beax max, Be. Then if we run Algorithm 5 with losses (SOR(()) - 1)/K normalized by K
log
Tog
for max = max, action set gli = 1 + (Wmax 1) for wmxx 2 1++
Breex
and context
1+V1-Pmms
discretization hyy Gmia C (j 1t). then the number of iterations will be bounded in expectation as
E
T
(WR)(W)) S mm + 4K 4K
( (CL /
1
log log at
+ [vr(w((a)) (29)
Waaxx
2
t=1
where L is the Lipschitz constant from Lemma E.I. Setting Whax
1
1t
Pmn
K
logg d
12g amas and m = 5/AC3L322T log" ami where & = 7 EE= 24CL
-Tog
log
and i. Amiin (D) + 4-min Amnia (D)+cmmu
149max
yields
Bwax
/1-B/
E
T
(w))
< / 54C3 amm + 4 log
1
24C17T3
+
(U,(w(()))
T
log² 2
ammax
1=1
(30)
5AC3133
4 log
s V log Otmax
2ACLTS
(i((((a)
log cimax
Proof. The bound follows from Theorem B.4 by noting that Lemma 2.1 implies that the functions
U, are -(1- log log" az (uaxx ) -semi-Lipschhizz over [1, Wmax] and the functions SOR; S U
are log o, -bounded. Note that for the choice of wmmx we the interval [1, wmax] contains the range of the optimal policy w", and further by Lemma E. it is L Lipschitz over [Cmin Guia + C].
E.4 EXTENSION PRECONDITIONED CG
While CG is an iterative algorithm, for simplicity we define as the solution to a minimization problem in the Krylov subspace:
Definition E.1. C,(A.b,w) Ain -blaSe le for X
arg min
|x 1--A-11 A
x=P. A/W
where the minimum taken over degree k polynomials Pk Raxn + Buxn
25

---

Lemma E.2. Let A be rasitive-definit matrix and b IR" any vector. Define
T llg
Vn(A)
x(A)
(w) = =+
(31)
-log
1 -
2+ I- p
Jus
for u Amax (DA- 2 1, Amax (LD- D/4)A-1) E [-1/4, 0), and the smallest constant (depending on A and b) s.1. U 2 CG(A,b, b, ). Then the following holds
1. 7 (0,1]
2. if Ps 1 then U is minimized at w" 1+ /1
and monotonically increases away
from wt in both directions
3. U
1-Aut1 Anv+2u-I TV TV ssmmi-iipschitz on [2v2- 2)
Flog( (((A)) on 12v2-22 1+1l T+I//Pmas |. where y S T+3pmmx
4. i S pmax then U
1+
log 1-1T5))
Proof. By Hackbusch (2016, Theorem 10.17) we have that the kth residual of SSOR SSORRprecondittonee CG satisfies
2x
||>(2))l22 =|lb- Axall VIAIIIA "b --xx/lA
||
A b
Xx A
1 222
(32)
2 K (<()
|roll2
+ 2
VE(W5 A)-1 for x
=
By Axelsson (1994, Theorem 7.17) we have
A)+1
VWW A)+1
1+ (2- ww
A)
(33)
2
Combining the two inequalities above yields the first result. For the second, we compute the derivative w.r.t. Wi
8(2v + 1)w² 4u((2 - w)²
(34)
7
(2 w)w
n(2 + (1122 w) 4w(ww w - 1))
2-
Since v E [-1/4,0] 0] and μ > 1. we have that p22 - w) 4w(w + w 1) 2 (2 - w) + 4w, which is nonnegative. Therefore the derivative only switches signs once, the zero of specified in the second result. The monotonic increase property follows by positivity of the numerator on wt The third property follows by noting that since y is increasing on w" only needs to consider [2v2- 2, W), where the numerator of the derivative negative; here we have
4p(2 w)
/ + 4u +
V 11V2
m(2-w) 4uw (n(2 w)2 4w(w - 1)) 41m + 2μ
(35)
where we have used (222- w) + Aw-v - 16w and w 2 2 - 2. For the last #+4u44
result we use the fact that K( (W.IA A) is maximal at the endpoints of the interval and evaluate it on +1)max+1]
those endpoints to bound < + (axa +,, /Pmas
8
We plot the bounds from Lemma E.2 in Figure 4. Note that E (0, 1] is an instance-dependent parameter that is defined to effectively scale down the function as much as possible while still being an upper bound on the cost; it thus allows us to exploit the shape of the upper bound without having it be too loose. This is useful since upper bounds for CG are known to be rather pessimistic, and are able to do this because our learning algorithms do not directly access the upper bound anyway. Empirically, we find often be around 3/4 larger.
26

---

E.4.1 PROOF OF THEOREM 2.3
Proof. By Lemma E.2 the functions U, - 1 2 CG2 - 1 are
M urv 2,22 -semi-
4ppm + 2p
Lipschitz and
log( /mix)
2 1744 -bounded on [2v2 - 2, 1+1/ note that by the assumption on
log Comma F
min; j4x and the fact that 1 2 1/4 the semi-Lipschitz constant is (l)))) Therefore the desired regret w.r.t any E [2v2- - 2, T+1/ /aaux follows, and extends to the rest of the interval because Lemma E.2.2 also implies all functions Ur are increasing away from this interval.
27

---

F SEMI-STOCHASTIC PROOFS
F.I REGULARITY OF THE CRITERION
Lemma E1. bla is p(C,,) ( ((A) (2-wmex + +2p(DA-)) 1) -Lipschitz w.r.t w E Q
Proof. Taking the derivative, have that
PuuliCablla () b
(² El o (ow Ow)
b
(8, (w)
(0w Cw 0
Σ 2 + -
< <1bbllv R())
= (
(36)
where the first inequality is due to third is due to the sub-muliiplicativity of the norm. The last line follows by symmetry of A Ab $ which implies that the spectral norm of any of power equals that power of its spectral radius, which by similarity is also the spectral radius of C... Next we use matrix calculus tool (Laue et al., 2018) to compute
the second is the triangle inequality, and the
1
2-
Og O, =
+
DDD/W L)- A(D + +L A(D/WWL')
w
A DDDD
2
ADD/#+L' -DD L)
w3
1 1 - -1 J2 AW D(D/W L)-
2 -- +
AW
(37)
so since
W 1 A 1/22 = -
1 C. < and
~
| (D/ + HIIDDW#LT)-- DW /
A
-WW A
I ||W D I. MalA
(38)
$ 2/ A A A
2p(DA I M. S 40(DA-1)
have by applying E [1, maax] that
|A A (a, A
2
+
2
-p(DA-1)
4
+ Sp(DA-1)
(39)
2
W
2 Wmax
n
28

---

Algorithm 8: Symmetric successive over-relaxation with an absolute convergence condition.
Input: E Raxn b E RRn, parameter w E (0, 2), initial vector x E R", tolerance e > 0 D L+ A // D is diagonal, is strictly lower triangular // compute third normal form // compute initial residual
WT + (D/w L)D (D/++(I)++L"))
To + b b-Ax
for k O.... do
if 1lrrll2
then return K
// return iteration count (for use in learning) X = W rk // solve two triangular systems and update vector Fk+1 b Ax
// compute the next residual
Output: k
Lemma F22 [[C[(c)bl) (c)bbl2 is AmialA)+cuul
w.r.t all c 2 Cmin
-Amin (c)),, where (c) denotes matrices derived from (c) = + cl
Proof. We take the derivative in the above proof of Lemma F.1:
10e (c)bllall = (() A (c) (clla)
(40)
We then again apply the matrix calculus tool of Laue et al. (2018) to get
2 - w (Do/
(c)
w A
w
+ w²
Ao/DD//>+417
w
(41)
2 Ac)DDC/ L)-
+
W2
2 (
A WW.
W
+ w
(0
w
By symmetry of (c) we have
A 1 (c [ (c)laa) IW [W
= IL, (L s s 2pp2221100))
(42)
Furthermore
|| (O - (c) M[()) (.)A (e)l2
(43)
S ||2 L. Mo(o)l S 4p(A-
and le 11W S 40(4 (c)
(44)
so by the lower bound of
A (c)) have the result.
Aunia Amiia
29

---

F.2 ANTI-CONCENTRATION
Lemma F.3 Let X E Ruxx be a nonzero matrix and = mu be a product of independent random variables m O and E I with m² E [0,n]] a x²-squared random variable with n degrees of freedom truncated to the interval [0, n| and u distributed uniformly on the surface of the unit sphere.
Then for any interval I = (5,E fore A> 0 we have that Pr[ E E
24 2.
Proof. Let f be the p.d.f. of b and be the p.d.f. of of'g ~ (N,., L.). Then by the law of total probability and the fact that b follows the distribution of g conditioned on | S we have that
/
Pr( I)
Pr[[ 4lb xx)
JI Pr
Xb|l2 l/b = xx)g(x)
J|x1i>n dg(x)
(45)
/
2
Pr([ Alg = x)dg(x)
1x1155n
< 2 Pr[[ E 1|g = x))g(x)= (lxxx)=2PPr((((() 2Pr( E I)
where the second inequality uses the fact that a x2 random variable with degrees of freedom has more than half of its mass below n. Defining the orthogonal diagonalization Q" AQ X X and noting that Qg ~ (O,,,, I,), we then have that
Xg =
Ap.ax 1=1
(46)
for i... XI...., ,n. ~ (O, 1). Let h, h1 and h_1 be the densities of Σ AR.IX An. and Eis AR.jX2 respectively, and let u(a) be the uniform measure on the interval (a, a + 2=A + A Then since the density of the sum of independent random variables is their convolution, we can apply Young's inequality obtain
Pr( Xglla 2 E I) = Pr (rg] Xgl3 E (E2, (E ()9))
20-2:4+22
< max a2=²
h(r)dr
max a2s?
u(x (-)(())(
||u * h b ([c,x))
||u * h1 h- I ((e,o))
< ||u * hn ((,,00)))) h- 1h-1llc
(47)
2=+122
< max a2c
In (x)dx
25
e ZApai
max
dx
o2c
2m 1p.4p
2(a + 2 A2)
2a
< max a2
A
A
xAA...]
A
Substituting into the first equation and using Ap. = |X|13 2 p (X) yields the result.
30

---

F.3 LIPSCHITZ EXPECTATION
Lemma F.4. Suppose b = mu, where m. and u are independent random variables with u distributed uniformly the surface of the unit sphere and m²² E [0, n) x²=squared random variable with degrees freedom Iruncated to the interval [0, n.. Define as in Corollary D.2, 3 = ming p(I, D Az), and SSOR(a) min (j k be the number of iterations convergence when the defect reduction matrix depends on some scalar x for some bounded interval X E. IfllC,b 32K 32K"Lx/2/5 -Lipschizz w.r.t. æ.
is LppCC, L-pCCh)---Lisschttt a.s. w.E.t. any æ E X then ESSOR is
34
Proof. First, note that by Hackbusch (2016, Theorem 6.26)
84
= =
= | A. 2 p(M)) 2
(48)
1+ /1 g2
16
30 /x/2
Now consider any I1, X s.t. | - -22 S
2K"L assume w.l.o.g. that #1 < 22, and pick
x' (2) with maximal p(Cz). Then setting Pz* pl ((zz) we have that ||0% bll2 is - Lipschitz for both i 1, and all E [K]. Therefore starting with Jensen's inequality we have that
[ESSOR(a;) ESSOR(2')
E[SsOR(I;) SSOR(x')
K
1 (r-SO(((SSO((:)) = &s sSSOR(2)) =
K
K
Cr bll2 en en //C E) + Σ
Pr(l bl12 s senl >
k=1
1>k
K
Pr(I E (E LP x =|)
K
K
Pr(1 Pr[ E (E,E Lp k-1 - - 2'll)
1=
1>
2 1-1 12/T p(O(.)
2 k- 12/
K
[x; 2'1+
x1 2'1
pl p(O²))
2K
LV 2/
32K L 2//T
Ij x'
Px
31
(49)
where the second inequality follows by the definition of SSOR, the third by Lipschitzness, and the fourth by the anti-concentration result of Lemma F.3. Since this holds for any nearby pairs æ x2. laking the summation over the interval N' completes the proof.
Corollary F.1. Under the assumptions of Lemma E.I. the function EbSSOR(A, b,w) is
32K (2nn(A)/=
2-wmaw 22((DA-1) -Lipschitz w.r.. w E [1, C (0, 2).
Proof. Apply Lemmas F.1 and E.I, noting that [lb||2 /ni by definition.
Corollary F.2. Under the assumptions of Lemma E.1, the function E, (,SSSR(A(C), b, w) is 320K' (2nx(((())) -Lipschitz w.r.t c 2 Cmin > 3=(Amin (A) +cmin)
max.
-Amin-
Proof. Apply Lemma F.2 and E.I, noting that [bll2 5 / by definition.
31

---

5000
Tsalis-INF
remainre
J= (mu')
Tsolis--NF-OO
instances remainig
11
Instanzes remainin minaiinn
11
ChebCB
Instance-Oplima
3000
6
Tsalls-NN
Tsalls INF
2000
2000
1000
3
2
6
10
12
14
0.5
1.5
2.5
3.5
1otal Berations
10tal Benations
5otal Tenations
Figure 5: Average across forty trials of the time needed to solve 5K diagonally shifted systems 12c-3 In for e Beta( 1 ?) (center) and Beta 6) (otherwise).
with A,
A+ 20
G EXPERIMENTAL DETAILS
All numerical results were generated in MATLAB laptop and can be re-generated by running the scripts available at https:/ /gi thub. com / mwhodak/ /aakk/leare metarninnto-ree since do not have access problem parameters, we experimented with few approaches to setting them automatically or heuristically on the simplest (low variance) setting below and then used the same settings for the rest of the experiments (high variance and heat equation). Furthermore, because the default step-size/learning rate settings in both algorithms are rather pessimistic, we use more aggressive time-varying approaches in practice. For Tsallis-INF we set 1ソt = 2/VL, which is an increasing
Note that,
what is used in the anytime variant (Zimmert & Seldin, 2021). As for ChebCB, schedule = (tt); note that Simchi-Levi & Xu (2021) also use an increasing learning rate schedule for setting inverse gap- gei-weighted probabilities.
G.1 BASIC EXPERIMENTS
For the experiments in Figure (center-left), we sample T 5K scalars e ~ Beta(2, (2, 6) and run Tsa,lis-NNF, Tsallis-INF-CD ChebCB, the instance optimal policy w (c), and five values of w evenlnly spaced on [1, 1.8] on all instances A, = A + 12c,-3 In, in random order. Note that Figure (left) contains results of the same setup, except with C Beta( 2. ). the higher-variance setting from Figure 1. The center and right figures contain results for the sub-optimal fixed w parameters, compared Tsallis
For both experiments the matrix A is again the 100 x 100 Laplacian of equare-shaped domain generated in MATLAB, and the targets are ra-saelled at each instance from the Gaussian truncated radially have norm S The reported results averaged of forty trials.
G.2 ACCELERATING A 2D HEAT EQUATION SOLVER
We then consider applying our methods the task of numerical simulation of the 2D heat equation
O,u(t, x) x(t) x) + f(t, x)
(50)
over the domain [0, 1]2 and E [0,5]. We use five-point finite difference discretization with size denoted n = 1/Axx so that when an implicit time-stepping method such as Crank-Nicolson is applied with timestep A, the numerical simulation requires sequentially solving a sequence of linear x((t + 1/2) for fixed matrix A (that depends on Ae
systems (Ac, b.) with A, I and Ax) corresponding to the discrete Laplacian of the system Veque, 2007, Equation 12.29). Each A, is positive definite, and moreover note that mathematically the setting is equivalent to an instantiation of the diagonal offset setting introduced in Section 2.4, since the linear system is A b, for c 1/m((1 + 1/2) However, for simplicity will simply pass x((t + 1/2A,)) as contexts to CB methods.
equivalent to e, I
To complete the problem specification, define the bump function (x, centered at c E g² with radius 0 be exp T- - if x c/22 < r and otherwise. We set the initial condition u(0, x) b 1). (x). forcing function f(t, x) 32b, +coo(166)))4 (oo.6r/))//./.o(11666/ (co0(166))4
and diffusion coefficient ro()) = max{0.01 sin sin())))), - 110 sin 2.tt). }. The forcing function effectivvel
bump circling around the center of the domain- chosen to ensure that the linear system solutions
32

---

vanilla CG
valdock CB waletokk
=F1
waldeokk CG walkdookk
1.5
=11
Tsallis-INF
diffusion coefficient
12
ChebCB
0.5
220
1253
01
10
10
10*
10
metit dimerrsionn
metrix dimerrsion
Diffusion coefficient as function of time (left) and normalized total wallelock time required to run 5K steps of the numerical simulation (center and right). The numbers within the middle plot corresponding to the average number of seconds required to run a step of the simulation using vanilla CG. The right-hand plot shows 95% confidence intervals across the three trials for Tsallis-NN and ChebCB at the three higher-dimensional evaluations.
Figure 6:
are not too close to each other or to zero, and the diffusion coefficient function- -plottedlteeddeee in Figure 6 (left)-is chosen to make the instance-optimal w behave roughly periodically (c.f. Figure 2 (right)).
We set A, 10-3, thus making T 5000, and evaluate approach across five spatial discretiza- tions: 7xx = 25, 50, 100, 200, 400. The resulting linear systems have size n = (nxx - - At each timestep, we solve each linear system using CG to relative precision consider are vanilla (unpreconditioned) CG and SSOR-CG with w = or w = 1.5; as comparators also evaluate performance when using the best fixed w in hindsight at each round, and when using the instance-optinal w at each round. Recall that we showed that Tsallis-INT has sublinear regret w.r.t. the surrogate cost of the best fixed w of SSOR-CG (Theorem 2.3), and that ChebCB has sublinear regret w.r.t. the instance-optimal w for SOR in the semi-stochastic setting of Section 3. Since both methods are randomized, we take the average of three runs.
10 S_ The baselines we
In Figure 2 (center-right) we show that both methods substantially outperform all three baselines, except 7 25 and 7 50, when w 1.5 almost recovers the best fixed parameter in hindsight; furthermore, ChebCB does better then the best fixed w in hindsight in most cases, In Figure 6 (right) we also show that-at high-enough dimensions- this reduction the number of iterations leads to an overall improvement in the runtime of the simulation. Several other pertinent notes include:
1. At lower dimensions the learning- based approaches have slower overall runtime because of overhead associated with learning; ChebCB in particular solves small constrained linear regression at each step. However, this overhead does not scale with matrix dimension, and we expect data-driven approaches have the greatest impact in higher dimensions.
2. Vanilla (unpreconditioned) CG is faster than SSOR preconditioned CG with despite having more iterations because each iteration is more costly.
3. To get a comparative sense of the scale of the improvement, we can consider the results Li et al. (2023, Table 1), who learn (deep-learning baaeddggg-aaadaeaeee preconditioner for CG to simulate the 2D heat equation. In the precision 10 case their solver takes 2.3 seconds, while Gauss-Seidel (i.e. SSOR with 1) takes 2.995 seconds, roughly 1.3x improvement. In our most closely comparable setting, INa and ChebCB are roughly 2.4x and 3.5x faster than Gauss-Seidel, respectively (and have other advantages such simplicity and being deployable an online fashion without pretraining). We caveat this comparison by noting that Li al. (2023) consider a statistical, not online, learning setup, and their matrix structure may be significantly different-it results from a finite element method rather than finite differences. The only way to achieve a direct comparisons is via access code; as of this writing it is not public.
Lastly, we give additional details for the plot in Figure (right), which shows the actions taken by the various algorithms for simulation n 100. For clarity all lines smoothed using moving average with window of 25, and for Tsallis-INF and ChebCB we also shade one standard deviation computed over this window. The plot shows that Tsallis converges to an action close to the best fixed w in hindsight, and that ChebCB fairly quickly follows the instance-optimal path, with the standard deviation of both decreasing over time.
33


Reviews:
[{"review_id": "XmxUsRbw7H", "reviewer": "Reviewer_9VfQ", "summary": "The paper considers an online parameter optimization problem for sequentially solving linear system equations with a specific iterative algorithm (SOR) with a parameter.  The goal is formulated as a regret minimization where each cost per trial is defined as the number of iterations with the parameter at the trial. The critical observation is that the cost is further upper bounded by a convex surrogate function and analyzed under the online convex optimization framework.", "questions": "Is it possible to extend this framework to other algorithms for solving linear equations?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem is well-motivated (even though the algorithm looks slightly restrictive). The critical observation is that the SOR iteration is bounded by a convex surrogate function with the parameter. The observation is non-trivial and thus the paper shows a new and interesting application of the online convex optimization framework.", "weaknesses": "Maybe a weakness of the paper is that the reduction to OCO is restricted to a certain type of algorithm (SOR) only so far. But I do not think the weakness is not so crucial since it shows, to the best of my knowledge, a new application of numerical optimization from OCO.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mZ2yt7sV82", "reviewer": "Reviewer_7bs2", "summary": "This a paper about learning a good choice of a parameter of a linear system solver in an online (bandit-like) fashion. The goal is to minimize the number of iterations (and, as a consequence, the running time) of the solver. The paper is theoretical in nature – the authors prove regret bounds; simple proof-of-principle-style experiments serve only a complementary role.\n\nThe solver of choice for the paper is Successive Over-Relaxation (SOR). It is not the most widely used in practice solver nowadays, but it has some uses, and it has the advantage of being simple and easy to analyse. The authors mention extending their theory to more complex solver as an important direction of future research.\n\nSOR is an iterative solver, and it has a single parameter, denoted by omega. The authors study the following setting: a sequence of linear systems is solved one by one. Before we start solving each system we can pick a value of omega, and after the system is solved we only learn for how many iterations the solver was running.\n\nThe authors propose to use Tsallis-INF bandit algorithm (a variant of more popular Exp3) to learn the parameter omega. The challenge lies in analysis – the number of iterations as a function of the parameter is non-Lipschitz. Instead, the authors use a continuous surrogate upper bound on the number of iterations, which they claim is reasonably tight. Unfortunately, this means they can only prove that the surrogate cost of their algorithm is close (up to a sublinear regret) to optimal surrogate cost (not the real cost, which might be much smaller) for a fixed parameter chosen in hindsight.\n\nThe second setting studied is in the spirit of contextual bandits. The authors assume all coefficient matrices are the same up to a linear shift of the diagonal. They show that under this assumption it is possible to obtain sublinear regret with respect to best choice of the parameter for each value of the shift separately. Since for a fixed shift there is a fixed optimal parameter, it means that the benchmark is just the optimal unconstrained choice of parameter. This is a very strong result (though only under a strong assumption about input, and still about the surrogate cost).\n\nFinally, the authors also study a setting where target vectors (but not coefficient matrices) are drawn from a distribution, and in this setting they are able to show sublinear regret bounds using the actual (not the surrogate) cost of the algorithm.", "questions": "Do you have any (even soft) argument that the limitations imposed (surrogate loss or stochastic input) are necessary to prove sublinear regret bounds?\n\nIt seems that learning optimal parameters for iterative methods (LP solvers, linear systems solvers, GD, etc.) should be similar and explainable with a single theory. Do you have any understanding why it does not seem to be the case, and we have separate results for each of these methods?\n\nHow your work compares to what people do in practice? Is it clear why they are not using similar methods? Would it be feasible to compare empirical performance of your approach to what is already done in practice?\n\nMinor remarks:\n\nAbstract: \"we prove that a bandit algorithm (...) can select parameters (...) such that the overall cost is almost as good as that the best fixed \\omega would have obtained\" – this sentence is technically not true, you prove it only for the surrogate upper bound and not for the actual cost\n\nPage 2: \"ir\" -> \"it\", \"in-addition\" -> \"in addition\"\n\nPage 4: \"known scalars c_t\" – I'm not sure what \"known\" means here; even if the scalars are not given explicitly, they can be trivially inferred from the input.\n\nPage 4: please remind the reader that rho denotes the spectral radius – it might be clear for optimization people but not necessarily for learning people.\n\nPage 14: \",e.g.\" -> \", e.g.\"", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The setting studied is very natural and this paper can easily stimulate further research in the area.\n\nThe algorithm used is a standard one, so the proposed approach seems more practical than if it was an ad-hoc algorithm designed specifically so that the analysis works.\n\nThe paper is nicely written – the authors explain why they do certain things instead of just presenting proofs out of the blue.\n\nIt seems that the technical content is novel and nontrivial – though I do not know the area well enough to be certain about that.", "weaknesses": "The results are either about the surrogate loss, or under the assumption that part of the input is stochastic and not adversarial.\n\nExperiments are very basic. It is not a huge weakness as this paper can be positioned as a pure theory paper. However, since the solver seems to be used in practice, it seems it would not be difficult to get some actual input data it is run on, and results of such experiments would be very interesting to see.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FFrjLxjCjU", "reviewer": "Reviewer_vVy3", "summary": "In this paper, the authors study different settings of sequentially solving a series of (related) linear systems using Successive-Over-Relaxation (SOR). The problem is that the runtime/number of iterations of SOR highly depends on the \"relaxation parameter\" $\\omega$, so the challenge is to learn online good choices of $\\omega$. They look at two settings. First, a deterministic setting where we assume the instances are close enough to the \"asymptotic regime\" where the convergence rate of SOR decreases linearly with respect to the spectral radius of a specific problem-dependent matrix. The second setting we assume the target vector is randomly sampled from what we can think as a truncated Gaussian distribution. In both cases, they show how bandit algorithms can be used to pick parameters $\\omega$ online that work well on average. Furthermore, they show how to use contextual bandits in the case when we have more structured information about the linear systems we need to solve.", "questions": "- A few times the authors mention that they provide \"end-to-end\" guarantees, but I am not sure if I can parse what you mean by this. Could your briefly mention that do you mean by this? \n- If the authors have time, could you expand a bit on the poor dependency on the spectral properties of the instances in theorem 2.3 (and maybe 2.2)? The authors mention that these dependencies are not ideal, but I do not have enough context to know what are dependencies that would be more \"acceptable\";\n- On the assumption on $b_t$'s in sec 3, is it really necessary for the entries to be independent of each other? Moreover, how much do you depend on $b_t$'s being exactly the distribution you have? I probably could answer this questions myself if I had the time to go carefully over the proofs, but if the authors could briefly comment on this, it would be great;\n\n### Suggestions\n- The discussion on \"asymptocity\" in sec 2.2 is very unclear. At some point I think I understood more or less what you meant, but this is a very confusing part early on in the paper, and there is not enough context to interpret the plot. Maybe this is more of a note, but if you could expand on this in a revised version of the paper, it would be great.\n- In Theorem 2.1 you cite \"Lemma 2.1.2\" to mean item 2 of Lemma 2.1, which is a bit weird since this is easily misunderstood as a separate lemma (that does not exist);\n- Although I mentioned about one of the plots in Fig 2 before, I think all of the plots in Fig 2 need more context. It is not clear to me what the middle plot is trying to show, and the last plot is great but I only managed to understand that all of the algorithms there were algorithms proposed in the paper by the end of the paper. So when I first read this plot I was very confused;", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "- This is a paper with a very solid contribution on data-driven algorithm design (and a few interesting technical bits for adversarial badits), with interesting new ideas applied to a simple but interesting problem, and I believe will likely lead to interesting follow-up work.\n\n    More specifically, generalizing bandits algorithms in the adversarial setting for functions that are not exactly Lipschitz continuous for the use in algorithms with predictions/data-driven algorithm design is very interesting and seems to work quite well. The extension to the use of Contextual Bandit algorithms when we are looking at a series of systems with shifted diagonals is interesting, and how we can exploit the stochastic case and skip the use of a surrogate function at all are quite interesting.\n    \n    Even if some of the bounds do not seem to lead to very informative regret guarantees, it seems to be full of interesting technical ideas and does a good job of showing what are the difficulties of controlling the dependency of problem-parameters in these regret bounds.\n    \n- As the authors themselves mention, this is an interesting case where they designed algorithms of two related lines of work (algorithms with predictions and data-driven algorithm design)\n\n- The authors are very careful to not over-sell their contributions: caveats and weaknesses are often immediately mentioned and discussed, together with a summary of these drawbacks at the conclusions sections;\n- The presentation is not perfect, but very good considering how much ground the authors cover in 9 pages. I believe this paper covered a lot of ground, and after skimming a few parts of the appendix it does seem to be very thorough\n- Although it is in the appendix due to space limitations, I really enjoyed reading the related work section since I lack the background in linear system solvers and the related literature;", "weaknesses": "- As mentioned by the authors, some of performance guarantees seem to be very loose (mainly the ones in section 3 with the dependency on dimension);\n- I am more acquainted with the online learning and optimization literature, and probably did not have enough background to judge how good some of the bounds were (mainly, the dependency on problem parameters in Thm 2.2 and 2.3). Although this is a problem with my lack of background, I believe the authors could add a bit of discussion of what would be \"ideal\" (?) dependencies here;\n- It seems to be that having stochastic target vectors with independent entries is a very strong assumption. The authors acknowledge that assuming the vectors come from a (scaled and somewhat truncated) Gaussian is restrictive, but it was not clear if mild correlation between the entries (positive definite correlation matrix with small off-diagonal entries?) would make the analysis break down;\n- Although the experiments are not the focus of the paper, I do believe they are a bit too simplistic. I do no think the authors should spend time coming up with more experiments, but I think that discussing the limitations of the experiments would be helpful.\n\n\nI am not an expert in algorithms with predictions/data-drive algorithm design (I'm more acquainted with online learning and optimization), so I will take into account the other reviews for my final scores. I might be willing to increase my score, the only reasons why I haven't said yet that this should be highlighted is because ICLR is a broad conference (with a slight focus on deep learning and other large-scale ML techniques), so I am not sure if I would say this should be highlighted, but the contributions are definitely substantive and solid.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "moWPbTL3pB", "reviewer": "Reviewer_655G", "summary": "This paper establishes learnability result for the overrelaxation parameter in SOR method. Online learning algorithm and regret bounds are analyzed for solving a sequence of linear systems.", "questions": "1. This paper focuses on online learning algorithms. In practice this may still be inefficient in the exploration phase. Is it possible to adopt a static learning model and prove learnability result using tools from statistical learning theory? (e.g., train a classifier/regressor to predict the best $\\omega$)\n2. The paper assumes iterative methods start from scratch. When solving a sequence of linear systems that are potentially from the same distribution, a common practice is to warm-start with previous solutions, and it often speeds up convergence. Is it possible also to incorporate this into your analysis?\n\n**Minor typos and stylistic issues**\n\n1. Page 14\n\n   guarantes => guarantees\n\n2. Page 26 \n\n   folloows => follows", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper theoretically justifies that algorithm parameter in linear system solving is learnable. Online algorithm with provable regret guarantees is given. The context setup and analysis are novel to my knowledge.", "weaknesses": "The paper adopts an upperbound surrogate loss function. Although this is reasonable due to the hardness in characterizing instance-specific convergence behavior of iterative methods, this still results in potential gap between theory and practice.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)