Variance Reduction of Stochastic Hypergradient Estimation by Mixed Fixed-Point Iteration
Authors: Naoyuki Terashita
Satoshi Hara
Venue: Accepted by TMLR ()

Abstract: Hypergradient represents how the hyperparameter of an optimization problem (or inner-problem) changes an outer-cost through the optimized inner-parameter, and it takes a crucial role in hyperparameter optimization, meta learning, and data influence estimation.
This paper studies hypergradient computation involving a stochastic inner-problem, a typical machine learning setting where the empirical loss is estimated by minibatches.
Stochastic hypergradient estimation requires estimating products of Jacobian matrices of the inner iteration.
Current methods struggle with large estimation variance because they depend on a specific sequence of Jacobian samples to estimate this product.
This paper overcomes this problem by \emph{mixing} two different stochastic hypergradient estimation methods that use distinct sequences of Jacobian samples.
Furthermore, we show that the proposed method enables almost sure convergence to the true hypergradient through the stochastic Krasnosel'ski\u{\i}-Mann iteration.
Theoretical analysis demonstrates that, compared to existing approaches, our method achieves lower asymptotic variance bounds while maintaining comparable computational complexity.
Empirical evaluations on synthetic and real-world tasks verify our theoretical results and superior variance reduction over existing methods.

Published in Transactions on Machine Learning Research (07/2025)
Variance Reduction of Stochastic Hypergradient Estimation by Mixed Fixed-Point Iteration
Naoyuki Terashita Hitachi, Ltd.
Osaka University
Satoshi Hara
satohara&@ucc.cu
University of Electre-Commuriiattinn
Reviewed on OpenReview: https: h/pp::///pperevveeee net/ formm?? mkmX2IC45e
Abstract
Hypergradient represents how the hyperparameter of optimization problem (or inner- problem) changes an outer-cost through the optimized inner-parameter, and it takes a cru- cial role in hyperparameter optimization, meta learning, and data influence estimation. This paper studics hypergradient computation involving a stochastic inner-problem, a typi- cal machine learning setting where the empirical loss estimated by minibatches. Stochastic hypergradient estimation requires estimating products of Jacobian matrices of the inner it- eration. Current methods struggle with large estimation variance because they depend specific sequence of Jacobian samples to estimate this product. This paper overcomes this problem by miring two different stochastic hypergradient estimation methods that use distinct sequences of Jacobian samples. Furthermore, we show that the proposed method enables almost sure convergence to the true hypergradient through the stochastic Kras- nosel skii-Mann iteration. Theoretical analysis demonstrates that, compared to existing approaches, our method achieves lower asymptotic variance bounds while maintaining com- parable computational complexity. Empirical evaluations on synthetic and real-world tasks verify our theoretical results and superior variance reduction over existing methods. The code is available al http::///ttthbb.occaaaaaaarrrree
1 Introduction
Hypergradient essential in wide range of tasks, including bilevel optimization which encompasses hy-
perparameter optimization (Maclaurin al. 2015
Franceschi et al. 2017 Liu al.] 2018 Lorraine al.
2019h and meta. learning (Andrychowiez al. 2016 Finn et al. 2017 Franceschi et al. 2018 well
explainability technique called data influence estimation (Koh & Liang 2017 Hara et al. 2019 Pruthi et al. 2020 Xue et al. 2021). Using a hyperparameter (or outer-parameter) À and an parameer æ that paramelerize a scalar cost f and vector mapping s, the hypergradient can be defined as
dxf(x(A),
with s<A) (((()),))
hypergradient
where dxf(x(A), A) represents the derivative of X (x)()),), which takes the differentiation of À + =(A) into account. A common choice for is to set p(x, A) x X) with 7 E R++ and a strongly-convex function g. In this case, =(A)= = ((x(() A), 3 A) is equivalent to the minimization problem s(A) = arg min, ((x, A).
Calculating exact hypergradicnts prohibitive for largo-scale problems as it requires the inverse of a parameter-sizer matrix, making their efficient approximations crucial. method called approximate im- plicit differentiation (AID) (rrazi et al. 2020 Pedregosa 2016 Lorraine al. 2019 Rajeswaran et. al.

---

Published in Transactions on Machine Learning Research (07/2025)
2019 considers the following form derived from the chain rule and implicit function theorem:
dxf( 02p(1 024) + 0xf.
where the arguments (a(A), are omitted and 0,f represents 0f/Oz. AID then replaces =(A) with the
inexact solution and approximately solves the resulting linear system. A common AID approach called fixed-point method [orraine et al. 2019) uses truncated Neumann series approximation with M terms solve the linear system,
(1- (_-))
Σ 0mip"
to avoid explicit calculation of the inverse matrix and Jacobiar
When the Jacobian matrix O2p is stochastically estimated, a typical scenario where the empirical loss is estimated by minibatches, computing hypergradient becomes even more challenging. To perform such an estimation, Koh & Liang] (2017) and Grazzi et al. (2021) used the stochastic fixed-point method which performs the following iteration with independent samples of 0zip denoted by (Am)mEN:
Wm+1 Am 0zf..
(1)
After M iterations, SM estimates the truncated Neumann series as:
M-1
M-1
M
= Σ
II
A
0of ~
Σ
m==
A=
m=0
This method can be seen as estimating the Jacobian product by,
Og 0. R AM- Am+ A..
for each m < M. This approach incurs large estimation variance because each order of the Jacobian product is sampled only once using specific sequence of Jacobian samples.
This work aims at reducing estimation variance by employing multiple different sequences of Jacobian samples to estimate any order of Jacobian product, without sacrificing computational complexity. We first point out that there are two hypergradient estimators that use diffcrent sequences of Jacobian samples to estimate their product: Stochastic fixed-point method (StocFP) (Koh Liang, 2017 Grazzi et al.] (022)) and Stochastic recurrent backpropagation (StocRB) (Ji al.) [021]] Yang al. (022)). Next, we propose MixedFP, which updates weighted average of StocFP and StocRB at each step. MixedFP realizes variance reduction by estimating the Jacobian product of given order using various sample sequences. Furthermore, we improve our MixedFP so it converges almost surely to the true hypergradient by applying the stochastic Krasnosel'skie Mann iteration (Grazzi et al. 2021 Bravo & Cominetti] 2024).
Additionally, this paper quantifics the variance reduction performance of the proposed method both theoret- ically and empirically. Our analysis showed that the proposed method improves the expected error compared to the existing stochastic bilevel optimization methods (Grazzi et al. 2021 Ji et al. 2021 Arbel & Mairal, 2022) in terms of one-shot hypergradient estimation. Numerical experiments using real-world data also sup- port the theoretical observations and demonstrate that proposed methods can estimate hypergradients more accurately than the existing methods.
The main contributions this work are as follows:
We point out that there are two algorithms in stochastic hypergradient estimation methods that estimate Jacobian products using different sequences of samples: Stochastic fixed-point method (StocFP) and Stochastic recurrent backpropagation (StocRB).
1 For simplicity of exposition, here we assume that the exact solution =(A) given. practice, ATD also approximates this solution.
2

---

Published in Transactions on Machine Learning Research (07/2025)
We propose MixedFP, which reduces the estimation variance by iteratively updating weighted average of StocFP and StocRB.
We achieve almost sure convergence to the true hypergradient by applying the stochastic Kras- Mkin-Maninn iteration to MixedFP.
Through theoretical analysis of error bounds, we demonstrate improvement over previous research (rrazzi et al. 2021 Ji et al. 2021 Arbel & Mairal 2022).
Our experiments using real-world data demonstrate a smaller estimation variance of the proposed methods compared to the existing approaches.
2 Related Work
Stochastic hypergradient estimation has been studied in the context of stochastic bilevel optimiza- tion (Ghadimi & Wang 2018 Coucllan & Wang 2016 Ji et al. 2021 to minimize estimation variances of solutions for both inner- and outer-problems.. Ji et al. 2021) h and Yang et al. (2021) propose methods that focus on two-time scale updates and stepsize adjustments to accelerate inner- pree--poolem optimization with a warm-up strategy. However, these studies do not address the variance reduction of stochastic hypergradi- ent estimation because the noise on the hypergradient estimation is managcable by decrcasing stepsizes of their outer-optimizationn. In contrast, our work specifically aims reducing the variance of hypergradient estimation itself. This is particularly important for tasks like influence estimation (Koh & Liang 2017; Khanna et al. 2019), which estimates the impact of excluding training data on performance by single-shot hypergradient estimation.
There are a few studies that addressed the variance reduction of hypergradient estimation. Grazzi et al. (2021) applied the stochastic Mosuoosl'sskiiMaane iteration (1), which is equivalent to solving min,, 1 1ll( d f by the stochastic gradient descent (Arbel & Mairal 2022), to achieve almost sure convergence to the true hypergradient. Moreover, Arbel & Mairal (2022) applies the warm-up to the hypergradient estimation to improve convergence properties. Our MixedFF incorporates their solution enabling a more general and effective framework for stochastic hypergradient estimation. Note that this paper does not restrict the context of hypergradients to bilevel optimization, and therefore the warm-up stratogy is not employed. Introducing a warm-up is a promising direction for future rescarch, which can enhance the convergence of the bilevel optimization solved by our MixedFP.
special case,
3 Hypergradient and Stochastic Approximate Implicit Differentiation
In this section, we first redefine the hypergradient and then introduce two existing methods for computing stochastic hypergradients. We also highlight the difference in the sequences of Jacobian samples used for estimation in these methods.
3.1 Hypergradient and Its Approximation
The hypergradient (*f() X) the total derivative of the composite function À IRdx + f(((A), A) E R, where f(x, A) E R is the outer-cost function and æ(A) e IRd is the inner-parameter solution. The solution of the tnmer-parameter z(A) is defined as the fixed point of mapping &p d x Rd + Rdd,, which gives the hypergradient as
dxf(x(A),A) with a(A) = ())
Approximate implicit differentiation (AID) Lorraine et al. 2019 Pedregosa 2016 is a set of approaches
consisting of the derivation of a linear system for (x)()(A), X) using the implicit function thcorem along with two approximations: an approximate inner solution x =()) and an approximate linear system solution. The implicit function theorem and chain rule rewrite the hypergradient in form that involves inverse
3

---

Published in Transactions on Machine Learning Research (07 (025)
matrix:
dxf(z(A), axp dzf 0 where d.f (T -T- -00
(2a)
(2b)
Here and hereafter, the arguments (z) are omitted when clear from the context. Moreover, for any vector function h 1 + R", denote its partial derivative by 0,h(2) E Rmxx, e.g,, 0xf represents the partial derivative of f(x(A), with respect to the second argument, treating z²A) as constant.
is justified under the following assumption. Assumption 1. For every Rdd,, we assume:
(i) A) is contruction; i.e., there erists constant < such that [10,p(a, A)ll S for any x e RR, (ii) p(x, A) and f(x, A) are differentiable at any x E IRd..
Since calculating the inverse matrix is expensive for large dz, Lorraine et al. (009)) proposed the fixed-point method, which approximates 2b using truncated Neumann scries:
d
Σ Oxpp"OO.. Oz
m==
This approximation becomes exact as M approaches infinity. Hereafter, we consider scenario where the exact inner solution z(A) is given, since our analysis in Appendix C Cindicates that the error from the incxact inner solution remains consistent among ours and previous work.
3.2 Stochastic Approximate Implicit Differentiation
From here on, we assume that we only have access to a stochastic estimator of ((x, A) denoted by 6(z, A; where is some random variable whose values lie within a measurable space For convenience, we introduce the following notations:
= ())) Á (=)
An As c= (=),)(x(AA))))
where (Em)MEN independent copies of &. We assume that & satisfies the following conditions: Assumption 2. For every x E gd. and X E Rds,
(i) ((x, E) is on unbiased eslimator of sps,, A); i.e. Elg(a, {: = p(a, A), Vx, A.
(ii) is differentiable with respect to the first and second arguments at any X.
Assumption i can be easily satisfied in typical machine learning scenario where the stochastic mapping
is a stochastic gradient descent, i.e., Ax = x - ~0,((x, A; and the instance & uniformly sampled from a finite datasct.
Note that, unlike common settings of stochastic bilevel optimization (e.g, Ghadimi & Wang 2018 ),
do not consider estimation of f. This is solely for clarity, as our primary focus is the estimation error caused by the inverse matrix approximation with stochastic & and our primary distinction from the previous studies Arbel & Mairal
2022
Grazzi et al. 2020 is the improvement, on this approximation.
3.2.1 Stochastic Fixed-Point Method (StocFP)
The studies by Grazzi et al. fixed-point method (StoclPP
2021 and Koh & Liang (2017) employ an iteration that we call Stochastic
stochastic variant of the fixed-point method (Lorraine et al. 2019
ii = c For m = O...., M 1:
(4a)
@m+1 Am
(4b)
4

---

Published in Transactions on Machine Learning Research (07 (075)
The iteration 4bb accumulates the products of Jacobian samples and their sum simultaneously. Assump- tion 2 guarantees that WMs is an unbiased estimator of the right-hand side of (3] A key advantage of this method is that the vector-Jacobian product Am E Id can be calculated in O(d2) time using reverse mode automatic differentiation [Baydin et al. 2018), because An Wm is the derivative the scalar function # + (0(,, X; (m) Wm evaluated al z(A). The total computalion time therefore O(Mdz) with a memory requirement of Odz)).
3.2.2
Stochastic Recurrent Backpropagation (StocRB)
Ji et al. 2021 independently proposed a method for estimating dzf using a different form iterations:
de = 0 R úo
(5a)
For m O.... M 1:
gm+1 Um + úm
(5b)
am+1 = Am
The final output ÚM is an estimator of d,f. Unlike StocFP 4bb,, this method divides the computation of (3)
into the product of Jacobians estimated by üm and their sum estimated by dm. We refer as Stochastic rocurrent backpropagation (StocRB) because this method
be understood backpropagation from f
with respect to æ that has recurrently passed through for M times.
While both StocFP 4 and StocRB (5) provide the same hypergradient approximation in expectation, they yield different estimates given the same samples. This is because they estimate the expected Jacobian product using different sequences of Jacobian samples:
Remark 1. For any m 1....,,, Am-1-... ... Ao to estimate the same value:
M, StocFP estimates A using ÂM- ... AM-m, while StocRB uses
Stochastic Fixed-Poin Point Method (StocFP)
Stochastic Recurrent Backprop. (StocRB)
WM = c + Am--1
(= A
jM =
(= A9c)
(< Atc)
Age
(~ A'c)
+ AAM-IAM-2
(N A?c)
A,Agc
(~ A?c)
+ AMMM1AM-2 Aie
(= (-1o))
AM-2 A1Agc (~ AM-1)
+ AMM-IAM2 ...
... Age AM)
+ AM-1AM-2 . A,Anc (~ A
4 Variance Reduction of Stochastic Approximate Implicit Differentiation
The proposed method is based on two ideas. First, in Section 4.1 we show that the estimation variance can be reduced by weighted average-like iteration between StocRB and StocFP. Second, Section 4.2 explains that almost sure convergence to the true hypergradient is achieved by using the stochastic Krasnoscl'she skii Mann iteration [Bravo Cominetti 2024). For clarity, proofs of all theorems and leminas are deferred the appendix.
4.1 Mixed Fixed-Pon Iteration (MixedFP)
Our MixedFP mixes StocFP and StocRB, motivated by the observation in Remark [ that the sequences of Jacobian samples used for estimation different from each other.
5

---

Published in Transactions on Machine Learning Research (07 (07/2025)
We introduce parameter called the mixing rate [0. 1] and derive the following MixedFP algorithm:
0p 0 Rd thg úg
(6a)
For m o,..., M 1:
0m+1 af6m ûmm) + (1 a))m
Wm+1 = Am + c
(6b)
um+1 dmdm
Here, im and ūm are the same as in StocFP 4 and StocRB respectively. 6m+1 estimates the hyper- gradient 3 and is equal in expectation to Wm and Im- The update iteration for 6m is slightly modified weighted average of tmm in (4b) and 9m in 45bp. In fact, one verify that 6m+1 tm when 0 and bm+1 Dm when a = 1. It is noteworthy that the update for 0m+1 is not simple weighted average of wm and Im. because uses the previously computed value im instead of fm- This recursive structure allows MixedFP to estimate Jacobian products using diverse sequences of Jacobian samples:
Remark 2. When a E 10,1) and given m E [I,...,MM, ŪM uses A ... Ax-m for any k such that m k k < M to estimate A More specifically, with some Ek am such that am,k > 0
Mixed Fixed-Point Iteration (MixedFP) 6
6 = c
(= Ad))
+
01,M-1 -1 AM-1 .... + 01,0 Ao c
(N A'c)
terms
+
02,M-1 AM- AM
+ 02.1 Å,
(~ A2c)
M terms
+
0M-1,M-
AM--1...AA A
-MM-1.M- AMs-2...
e
(2 (2 1M-1
+ AM-1AM-2....
Ag
(MAM))
This diversity in the products of Jacobian samples cannot be obtained by simply taking a weighted average of im and im. While MixedFP maintains the same order of computational complexity as StocFP and StocRB, it requires two nackor-Jacobian-proddee
(Ám um and âm each step. This results a higher per- iteration computational cost compared to StocEP and StocRB, which only requires single JVP computation. However, our experimental results on computational cost. in Appendix D.3 show that MixedFP remains competitive with existing methods even when evaluated on a wall-clock basis.
Under the following regularity conditions, we show that the expected error of MixedEP improves over StocFP and StocRB.
Assumption 3. For every À E Reda, we assume that:
(i) (II < for every E gd.
(ii) There erists a constant P < such that [|0, x: S P for every x. E TRd and & E X.
(ii) The function I(,2)) is Lápschitz continnous with constant I, 2 0.

---

Published in Transactions on Machine Learning Research (07/2025)
Method
Base estimator
KM iteration
Expected lo error bound
Grazzi al. Ji et al. [2021
2021 with
StocFP
01 + O(PT))
StocRB
01 + Opp")
Fi
StocFP & StocRB
(14 [4 ) 101++0 Opp"))
1 1+00
Grazzi al. 2021 with constant
StocFP
on + O(-)))
MixcdFP-KM with constant ") StocFP & StocRB
1+0. (I+va,) on + (pp)))
Table 1: Comparison of the upper bound of the expected & error in stochastic hypergradient estimation.
Our MixedFP and MixedFP-KM can reduce non-decaying errors by [Itva,) < 1 times smaller than their
140g
countcrparts. All symbols are defined in Thcorems 4.1 and and we assumc a S ý this table.
Theorem 4.1 (MixedFP). Suppose Assumptions to hold and [0, then for any m 0,
σ1 + (pf")
if E (0,1) otherwise,
E [110m dzfl1 -
1 1+a1
+ +0 a1)"")
where
212 01 = (1- (1-) PI p==p² an a
1+a This result shows that the scale of the non-decaying term, (I S 1, becomes smaller as a increases, but when or is too large, more specifically when > A. the decay rate becomes slower than O(p"). Recalling that a E {0, }} corresponds to StocFP and StocRB, this indicates that chosen in < < ý improves the non- decaying term and the decaying term over these conventional methods. Table D summarizes these findings. This analysis also suggests that for sufficietlly large m, the smallest error is achieved when a is infinitely close to, but strictly less than, one. However, for finite m, smaller a may be optimal if accelerating the decaying term outweighs the reduction of the non-decaying term. These observations align with the empirical results in Section 5.1
4.2 MixedFP-KM: Application of the Stochastic skai-nnn (KM) Iteration
This section shows that the stochastic KM iteration can applied to MixedFP, enabling almost sure convergence to the true hypergradient.
The stochastic KM iteration (Grazzi et al. 2021 Bravo & Cominctti 2024 is an algorithm that finds the fixed point of contraction mapping using its unbiased estimation:
Theorem 4.2 (Stochastic KM iteration (Grazzi et al, 2021 ). Let Ç be a random variable with values in a measurable space Z and (Sm)meN be independent copies of Ç. Let T Rd - Rdd be contraction mapping and suppose a mapping f. R2d x - Rdd satisfies E[T (s: [)]= T(s) and v[;(s; s si + s2) a|| for any Rd. Consider the stochostic Krasnosel skił- Mann iteration defined by
8m+ (1 7m 7mm)mm )m T (m))
(7)
where stepsizes (Mm)mEN satisfy
m
VKEN,
Ak
Σ Th
J0
Σ m
00.
(8)
I+ 52
k=U
Then the sequence (mm) converges almost surely to the fixed point of T.
7

---

Published in Transactions on Machine Learning Research (07/2025)
Next, we show that MixedFP is contraction mapping. To do this, reformat MixedFP into the following mapping of concatenated vector.
(x; Bz + d.
(9)
al al
(1 (1-00 a
0 d= E 3dz
where
w E Rddee
B = 0 A
0
E R dde x3dz
24
0 Å
0
The mapping F has the following favorable properties that enable the application of the stochastie KM iteration.
Lemma 4.3. When Assumptions and 2 mapping that has the unique fixed point [d.f" def" 0 IR 3d2.
hold, the function F(z) ([z))=B(((2:)) [[f(2;)) with a E [0, 1) is a contraction
By applying 7 to the mapping F in 9 we obtain the following algorithm, which we name MixodFP-KI
vo = Rd. Wp = uo =
For ..... M 1:
Um++1 = (1 1/m) Um + Tmm (a(Um + um) + (1 - a) a))m)
Wm++1 (1 7m Wm + 7mm (m + c)
7mm + Mm Am Um
Um++1
From Theorem 4.2 and Lemma 4.3 it follows that um converges to d with appropriate scheduling of Tmm-
Theorem 4.4. Let 1), max pt, and q maxfa, P}. When Assumptions
3 hold and mm
satisfies 8 with 22 2(²² + /²1 -)², then
lim Um = dzl a.s. m=+00
For fixed stepsizes, the error can be bounded by an exponentially decaying term and non-decaying term. Theorem 4.5 (MixedFP-KM). Let a E 1), When 7m < T-p and Assumptions to hold, then for any m 2 0,
1+a an
E [l - dcSl113 <
nn + O (t)
(1 van)
where
21y13,0²
on (1 p)- 711 p1) p)¹
Pn (1 np) on (1 + ma)3.
For comparison, we show our result for the existing method (Grazzi et al. 2021 which applies the stochastic KM itcration to StocFF Theorem 4.6 (Stochastie KM iteration of StocFP). Suppose 1mn = I and Assumptions to 3 hold.
Then for any 0,
E [[ - d < o, (pT))
where on and on are as defined in Theorem 4.5
The non-decaying terms in both Theorems 4.5 and 4.6 can be made arbitrarily small by the choice of v) at the cost of slow decay rate. The observations obtained from comparing these two theorems are essentially the same as those obtained in Theorem 4.1 That is, when 0 < a < é. both the non-decaying and decaying term This was also empirically confirmed in Section 5.2
improve over the existing method (Grzzi et al. 2021
We note that our theoretical analysis for decreasing stepsizes (Theorem 4.4 provides only asymptotic conver- gence guarantees without non-asymptotic error bounds. While our experiments employ decreasing stepsizes, demonstrating its practical performance, the fnito-time performance analysis is limited the fixed stepsize case (Theorem 4.5). Future work could focus on extending our theoretical analysis to provide non-asymplolic error bounds for decreasing stepsizes.
2Theorem 4.6 differs from result obtained from Grazzi solely for fair comparison with method, and some cases
2021
Theorem 4.2). discussed in Appendix B.2 this is
tighter than the original bound.

---

Published in Transactions on Machine Learning Research (07 (0/5)
y= 10-3
y= 10-2
y= 10-1
V=
error from Vxf
100
500
1000
0
500
1000
0
500
1000
500
1000
Number of iterations
Number of iterations
Number of iterations
Number iterations
h=0 (StocFP)
a=
a=0.99
a= 0.999
a=1 (StocRB)
Figure Mean L2 error of hypergradients estimated by MixedFF using different values of a. We adopted synthetic inner-problem p(z(A), -I-- (A) BA and plotted results obtained with different 7
5
Experiments
In this section, we first investigate the relationship between the newly introduced parameter a and the convergence properties of MixedFP. Then, we compare the hypergradicnt estimation accuracy of our proposed method with existing approaches various real-world task settings.
5.1
Effect of Mixing Rate
5.1.1
Settings
This experiment evaluates the hypergradient estimation error in a synthetic setting where explicitly control the Jacobian matrix samples. Specifically, we compute the hypergradient for the following case:
f(x(A). = = x(A) + dd= x, p(x, A) = (I - BX,
where A E Rdaxds is a random variable sampled from discrete uniform distribution Á H.)). Each matrix in HI...., H, was constructed such that its eigenvalues follow the uniform distribution over [0,1- c| with a small constant e E R++ to meet our assumptions. We tested differ- ent y over (0, 1], which corresponds varying the values of p and ý. The other coefficients E RRd,, d E RRd,, and B xdd were generated by sampling their elements independently from uniform distribution over 1). Note that and vanish upon differentiation, thus their values do not affect the hypergradient estimation. The experiment run 100 times with different seed values for sampling H.
5.1.2 Results and Discussion
Fig. 1 shows the mcan squared error of hypergradients at different y and a.
0.999 consistently outperform both a 0 (StocFP) and = (StocRB) across different values of 7, highlighting the benefit of our mixing strategy. The results suggest that the optimal value of a is close to but not exactly one, aligning with the observation from Theorem 4.1 The superior performance a = 0.99 over the larger a = 0.999 is also consistent with Theorem 4.1 for finite iterations, a slightly smaller can be optimal when faster decay outweighs the increased non-decaying term.
Estimations with a 0.99
Additionally, we observe that increasing 7 accelerates the decay rate while enlarging the non-decay term. Higher values of y result in smaller values for both p and á which aligns with theoretical finding in Theorem that implics smaller á leads to larger non-decaying term and faster decay rate P1.
Furthermore, some settings show interesting patterns in their error curves. When 7= 10-1 with a 0.999 or 0.99, the error seems to decay at two distinct rates: initially faster, then slower after a certain number of iterations (around m = 50). This seems reflect the bound 0(maxfpi,ar )m)) which yielded from the sum of two terms with different decay rates P1 and aj (as in Appendix A.5 39 the term with the faster decay rate converges quickly, and for sufficiently large m, the slower decay rate eventually dominates.
9

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Synthetic Problem
10a
StocFP
12 error from V
StocRB
MixedFP
10
StocFP+KM- (StocFP KM-iteration with constant n)
StocFP+KM1 (StocFP KM-iteration with decreasing n)
100
MiXedFP*KM- (MixedFP-KM with constant n)
MixedFP+KM (MixedFP-KM with decreasing
0k
2k
5k
7k
10k
Hyperparam. Optimization
Influence Estimation
Meta Learning
43 error from V,F
10 1
10-2
10
10
10-3
10-3
10k
0k
2k
5k
7k
10k
2k 5k
7k
10k
Number of iterations
Number of iterations
Number of iterations
Figure Comparison of stochastic hypergradient cstimation methods evaluatod on different tasks. Each line shows mean squared error over 10 trials with different Jacobian sample sequences (Am) Am.me Method- specific parameters were optimized by grid search, selecting settings with minimal average over the last 1,000 iterations.
5.2 Comparison with Existing Approaches
In this section, we consider several real-world machine learning settings compare the hypergradient estima- tion accuracy of our proposed method with existing approaches across various combinations of the dataset, inner-problem, and outer-cost.
5.2.1 Tasks
We evaluated the performance of hypergradient estimation methods for the following task settings:
Synthetic problem evaluates the estimation accuracy of hypergradients of the synthetic outer-cost with the synthetic inner-problem used Scction 5.1
Hyperparameter optimization performs binary classification on Adult Income dataset (Becker
Kohavi 1996), where the outer-parameter is a vector of regularization coefficients, each corre-
sponding dimension of the inner- parampparetter.
Influence estimation conducts multi-class classification on Fashion-MNIST Xiao et al. 2017) where the unter-paramett is a vector of loss masks, each assigned to a corresponding sample.
Meta lcarning tackles regression problem on the California Housing datasct PPace & Barry 1997), where the outer-parameler delermines a biased regularization term.
Except for the synthetic setting, we adopt Φ(x, λ; = x000g X; E) and define f and for each setting as follows. We denote training inputs and outputs by Ei and Eout, and use ( E E for the independent validation dataset.
Hyperparameter optimization refers to the hypergradient estimation problem in a bilevel problem where the regularization coefficients hyperparameters. Specifically, in Grazzi et al. (2021), we compute
10

---

Published in Transactions on Machine Learning Research (07/2025)
the hypergradient for outer-probles that finds the optimal hyperparameters À whose elements different regularization coefficients for each dimension of x:
f(z,1)=
Σ (Ein Soon
BCE(S Sin 2),
g(z, A: Coun)) = BCCE(fout <m + 2² diag(A)x, (10)
where Lin E RRM, Sout E {0, 1), BCE(p, g) = log (1 + (xp(-pa)), and diag(A) denotes the diagonal matrix whose diagonal entries are the elements of X.
Influence estimation (Koh & Liang 2017 Khanna et al. 2019) estimates the change in performance when sample is excluded from the inner-problem. Influence estimation is equivalent to estimating the hypergradient of outer-parane petee-parameter that determines the loss mask of each sample. That is, with the sample index j = Uniform((1....; N)) with dataset size N,
f(z 1)
Σ
ClECout CH x),
((1)):1) A, A,C.EEouut Em)z),
(11)
(E )
where CE(p, q) - Pe log (cxp(g))/Et) (xppp=)))== and (Lim.j Sout, R {0, 1) denotes the j-th sample in the training dataset. For the values of the loss mask, let E {1} be the ones vector and 1fj 10, denote vector with ones except for zero the j-th index. Then, we can express the inner-solution obtained including excluding the j-th sample sx1) or x(14)), respectively. The linear approximation of the change in f when the j-th sample is excluded from the training dataset is given by
f(x(1+)), 144) - f(x(1),1) ss ())( (14) - 1) = dc, E..
By computing a vector - dx.( E RN, we can obtain this approximation for all samples in the dataset. Thus, influence estimation reduces to estimation of the hypergradient.
Further details on influence estimation and its connections to existing work (Koh Liang 2017 are provided
in Appendix D.2.2
Meta learning is the following task adopted in Grazzi et al. 2020), which simplifies the task called implicit Rajeswaran et al. 2019).
meta lcarning Denevi et al. 2019
(z,
Σ
x(A) - Eone 11
g(z, A; (Eins Cout)) 21)))2l1e -1 Eour 11 * 11x l/c All?, (12)
1
(Ein Foun) EE
where (Sin, (ont) E 8 x XR PRE and A E IRd..
5.2.2 Baselines and Procedure
Hypergradient cstimation is performed using the five methods: StocFP, StocRB, MixedFP, StocFP with KM-iteration (denoted StocFP+KM+ and StocFP*KM4), and MixedFP with KM-iteration (denoted MixedFP+KM+ and MixedFP+XMJ). KM+ and KK4 represent constant and linearly decreasing 7m. respectively, and we used Tm (6 + m) as the decreasing stepsize following Grazzi et al. (2021).
For each setting and method, we conducted the experiment 10 times using different seed values to vary the sequence M, and calculated the mean squared error & with respect to dzf. We computed the exact values of 0zʃ, Oxf, and 0xy in 2 ensure that the only source of randomness came from (Am)m<s M We approximated z(A) as accurately possible by minimizing g using full-batch optimizations until they converge. The impact of approximation errors in z(A) is analyzed only theoretically in Appendix C For the outer-parameter X, we used task-specific initialization values without performing onter-optimization.
To ensure a fair comparison between the baseline and proposed methods, we use the parameters that yield the best performance for each method. These parameters include 7,, a, 3 and 8, which were determined vigrida search by minimizing the average crror over the last 1,000 steps.
More detailed settings deferred to our appendix.
11

---

Published in Transactions on Machine Learning Research (07 (025)
5.2.3 Results and Discussion
In all tasks, MiXedFP+KM4 yielded the best performance (Table D. Our theoretical analysis quantifies the improvements achieved under a fixed stepsize (by comparing Theorem 4.5 and Theorem (.6), while for a de- creasing stepsize, we have only established asymptotic convergence (Theorem 4.40. Nevertheless, these results demonstrate that, in practice, MixedKM can also improve ostimation accuracy with decrcasing stepsizes. When using KM-iteration with constant slepsizes, MixedFP+KM- achieves lower error than its counterpart StocFP-XM-, which supports our results of Theorem 4.5 Among methods without KM-iteration, StocFP and StocRB tend to reach plateaus in the early iterations, while MixedFP continues to reduce error over a relatively larger number of steps.
6 Conclusion
We presented a stochastic hypergradient estimation method that reduces variance incurred by estimation error of the product of the inner-iteration's Jacobian matrices. Our proposed MixedFP combines existing StocFP and StocRB estimators to leverage their distinct sequences of Jacobian samples, thereby reducing estimation variance without increasing compulational complexity. Moreover, by applying the stochastic KM iteration into MixedFP, we established almost sure convergence to the true hypergradient. Empirical evalu- ations on synthetic and real-world tasks, including hyperparamcter optimization, data influence estimation, and meta learning, validated our theoretical findings and demonstrated improved estimation accuracy in practical scenarios.
Future work will focus on extending our analysis to the case of docrcasing stepsizes and exploring efficiency improvements in the application, including bilevel optimization.
Acknowledgements
This work was supported in part by Grant-in-Aid for Scientific Research(B) 23K28146.
References
Marcin Audrychowicz, Misha Denil, Sergio Gomez, Matthew W Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando De Freitas. Learning to learn by gradient descent by gradient descent. In Advances Neural Information Processing Systems, pp. 3981-3989, 2016.
Michael Arbel and Julien Mairal. Amortized implicit differentiation for stochastic bilevel optimization. In International Conference on Learning Representations, 2022.
Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey A. Radul, and Jelfrey Mark Siskind. Aulomatic differ- entiation in machine learning: a survey. Journal of Machine Learning Research, 18(15)):1-43, 2018.
Barry Becker and Ronny Kohavi. Adult. UCT Machine Learning Repository, 1996. DOI: 10.24432/C5X3 C20
Mario Bravo and Roberto Cominctti. Stochastic fixed-point iterations for nonexpansive maps: Convergence and error bounds. STAM Journal on Coutrol and Oplimization, 62(1):19112)) 219, 2024.
Nicolas Couellan and Wenjuan Wang. On the convergence of stochastic bi-level gradient methods. Opti- mization, pp. 13833, 2016.
Giacomo Denevi, Carlo Ciliberto, Riccardo Grazzi, and Massimiliano Pontil. Learning-to-learn stochastic gradient descent with biased regularization. arXie preprint orXiv:1903. 10399, 2019.
Chelsea Finn, Pieter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In International Conference on Machine Learning, pp. 1126-11231 1125,113,,111 2017.
Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil. Forward and reverse gradient- based hyperparameter optimization. In International Conference on Machine Learning, pp. 1165-1173, 2017.
12

---

Published in Transactions on Machine Learning Research (07/2025)
Luca Franceschi, Paolo Frascomi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil. Bilevel pro- gramming for hyperparameter optimization and meta-learning. In International Conference on Machine Learning, pp. 1563 1572, 2018.
Saeed Ghadimi and Mengdi Wang. Approximation methods for bilevel programming. ar preprint arXiv:1802.0222: 2018.
Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. arXie preprint arXiv:2006.1622 2020.
Riccardo Grazzi, Massimiliano Pontil, and Saverio Salzo. Convergence properties of stochastic hypergradi- ents. In International Conference on Artificial Intelligence and Statistics, pp. 3826-3834 2021.
Satoshi Hara, Atsushi Nitanda, and Takanori Maehara. Data cleansing for models trained with sgd. In Advances in Newral Information Processing Systems, 2019.
Kaiyi Ji, Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning, pp. 4882-4892, 2021.
Rajiv Khanna, Been Kim, Joydeep Ghosh, and Sanmi Koyejo. Interpreting black box predictions using fisher kernels. In International Conference on Artificial Intelligence and Statistics, pp. 3382-33999 PMLR, 2019.
Pang Wei Koh and Percy Liang. Understanding black box predictions via influence functions. In International Conference on Machine Learning, pp. 1885-1894, 2017.
Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: Differentiable architecture search. arXiv preprint arXiv:1806.09055 2018.
Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit differentiation. arXiv preprint Xix::1911. 02590, 2019.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Cradient-based hyperparameter optimization through reversible learning. In International Conference on Machine Learning, pp. 2113-2122, 2015.
R Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(: 291 297, 1997.
Fabian Pedregosa. Hyperparameter optimization with approximate gradient. In International Conference on Machine Learning, pp. 737 746, 2016.
Garima Pruthi, Frederick Liu, Satyen Kale, and Mukund Sundararajan. Estimating training data influence by tracing gradient descent. In Advances in Neural Information Processing Systems, pp. 19920- 199330, 2020.
Aravind Rajeswaran, Chelsea. Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. In Advances in Newral Information Processing Systems, pp. 113-124, 2019.
Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv: 1:Xiv:170.00771 077477 2017.
Yihao Xue, Chaoyue Niu, Zhenzhe Zheng, Shaojie Tang, Chengfei Lyu, Fan Wu, and Guihai Chen. Toward understanding the influence of individual clients in federated learning. Tn the A AAAT Conference on Artificial Intelligence, pp. 10560 10567, 2021.
Junjie Yang, Kaiyi Ji, and Yingbin Liang. Provably faster algorithms for bilevel optinization. In Advances in Neural Information Processing Systems, pp. 13670-13682, 2021.
13

---

Published in Transactions on Machine Learning Research (07 (025)
A Theoretical Analysis of MixedFP-KM
In this section, provide detailed analysis of the proposed MixedFP-KM. We first recall the iteration schemes, then present key lemmas giving explicit forms of the iterates and the corresponding error bounds.
A.1 Recap of the Iterations
Let 77 E R++ be stepsize. We denote by (Um), (Wm); and (um) the sequences generated by the following iteration:
to 0,
(13)
W/ =
(14)
Ug =
(15)
and for m = 0, 0,,,
mm+1 (1 11 77 (a (um + Um) + am p' + 17 (cetmm + @W))
(16)
Wm+11 (1 11) + 1) (AmWm +
Pm + 17c,
(17)
Um+1 (1 17) + 17 (AmuU)) Pm um
(18)
where a E [0. a= a, r= == 1-=------ 17 (1 - a), and Pm -I- Am)
A.2 Proof of Asymptotic Convergence of MixedFP-KM
This section provides the proofs for Lemma 4.3 and Theorem
Lemma 4.3. When Assumplions and 2 hold, the function F(z) ([F(== with a E [0,1) 1) is contraction mapping that has the unique fixed point d def 0- 3d.
Proof. Let B = B.B). Because B is block triangle matrix, its eigenvalues are those of diagonal blocks, that are ensured to lie in [0, 1) by Assumption unique fixed point.
1 and E [0,1). 1). Hence, F is contraction mapping which has a
Let the first, second, and third blocks of the fixed point = = T(²") be v", w", and u". From Assumption we know that u² = Under Assumption 2 w² satisfies w" Aw + thus w² = (1- A)- --- da.. Finally,
since we know that v* av (1 f. which is true only when v* d&f since 1 a #0 0.
Theorem 4.4. Let a E 1), max(ax (a, p}, and q max(a, P}. When Assumplions 7 hold and ym
satisfies 8 with 52 = 2( + 42) - (),, then
lim Um = dzf a.s. m=
Proof. Let = [B]] and q = || where B = E Since is a block triangular matrix, its cigenvalues are determined by those of its diagonal blocks, giving us g = max p} from Assumption [1] Similarly, we
can establish that g max P) under Assumption 3. Applying Theorem 4.3 from Grazzi et. al. (2021).
oblain 2 Therefore, combining Theorem 4.2 and Lemnna 4.3 conclude that the statement (I-q]
holds true.
A.3 Closed-Forre Expressions of MixedFP-KM
We first derive closed-form expressions (or general forms) for the iterates (Um), (wm); and (um) under the MixedFP-KI iteration [16) to (18) with a fixed step-size 7). Note that in this section, we exclusively consider the case where the step-size remains constant throughout the iterations to avoid any confusion.
14

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Lemma A.1 (General-Toee Term Representation). Under the MizedFP-KM. iteration slep-size mm, for 2 0 we have:
13
to
18
with fized
Um++1
m- (aau + cork),
(19)
R=0
Wm++1
I
C
Ps
(20)
1=0
k=0
j=++1
Um++1
II
Pj
C
(21)
1=
where E [0,1], ā == a, -
77 (1- a), a), and Pm 1 17 Am)
Proof.
We prove each identity by induction.
(Proof of 19 For m = 0, 166 matches
19 since
01 1/r (aro + awo)
Assune
19
holds for some m 2 0. Then from 16
U Um+2 = r1m+1 + (oul+11 + am
r m (orux + awns) + 17 (aum+1 + (10++1)
r
k=
m++ (orua + oww)).
k=U
This completes the inductive step.
(Proof of 20 For m 0, from
17 have
W1 Poe + njc,
since we = c. This is consistent with
20 for m = 0 because
W1
II P
c
c+1
II j=k+1
PS c Poc + nJc.
1=0
k=
Assune 200 holds for some 2 0. Then from
17
Wm42 Pm+ Wm++1 17c
m
m
772
Pm
II
P3
c
P
Pm+t1
Σ
II
Ps
c ne
A=U
j
m++1
II
P
c+
P3
j=k+1
This completes the proof for Wm+11
(Proof of 21
The derivation is analogous to that of Wm, except that no extra term added. For 0. 21 and 88 are consistent because both are
W1 Poo.
15

---

Published in Transactions on Machine Learning Research (07 (07/202) /025)
from Mo
The inductive step is essentially the same as above, giving
m++
Um+2 P
um++1
II j=0
Pj c.
Hence
21 holds.
A.4
General Error Terms
This section considers differences wetween (Um; Wm: Um) and their deterministic counterparts (Um; Wm: Um)) namely,
Uo = 0. To = C,
(22)
(23)
To = C,
(24)
and for m = 0,,,2,
mm+1 (1)1)) 11) Um (ax (Tm 7m) (m) mu (0 mm))) Um+1 (1 17) Wm 17 (Awn m + c) Pm Wn 17C, um = (1 - 17) + (Aum) Pm Um
(25)
(26)
(27)
where P= () A). The next lemma gives a useful expression for these differences.
Lemma A.2 (General Error Terms). Under the MixedFP- PMM iteration 13 to 18
with Tm = 77 and its
deterministic counterpart 22 to 27 for s 2 1 we have,
8-1
U. 0 1
(a (uk - Tu) + + (we wh))
(28)
ps-1- (Ag -A) {
W.
W 1)
II
P
c +
Σ
II
I
(29)
1=
1=0
k=0
Ni=k+1
8-1
k-
U1 T.
ps-1-k
(Au - -A)
II P
(30)
1=0
where P nyT - A).
Proof. We give the argument for each the three differences, using induction and the results in Lemma A.1
(Difference for Um) As the direct consequence of Lemma A.1 the general term of the deterministic case is given by
Uo+1
=1]
(aTT + awk),
(31)
k=0
One can directly obtain 28 by comparing the general terms 19 and 31
(Difference for w.) By subtracting (26
from
17 the following recursion is shown to hold:
Ws+1
Wst1
P.w. Pw, + n/c n/c
P (w,, W.) + (P: P)ws
P P(w. W.) n(A, A)wa.
(32)
16

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Hence, 29 is consistent with 32
when 0, because both are
W1 1 -n(A0- A)c,
(33)
using the fact that wi = Wig = c. Assume 29 holds for some & 2 1. Then, from
32
and
20
Wa+1 1s+1
P (w,- w. + [1A A)wa
s-1
1-1 (A, A) {
5-1
np _ po-1-1
III, P
P
0
+ +1 1)
II P
1=0
k=0
pa-j (A, - A)
II
P
c+1
1
j=0
k=
Li=k+1
This completes the proof for +
(Difference for u.) Similar steps yield
30 without the extra term in the iteration for um+11 Namely,
substraction of (18) from [27 yields
U To (u, T.) + n((A A)ua,
and one can show that this consistent with 1, then 30 holds true for + because
30 when - O. When 30
assumed to be true for some
Us411 T941 = P (u, - T.,) + A(A A)us,
(34)
np
ps-1-i(A, - -A)
II
P:
c ( - A)
II P,
(35)
i=0
i=
j-1
p
(A, A)
II
P
(36)
1=0
i=0
A.5 Bounding Variance Norm of MixedFP-KM
We prove the bound on E [|um [ starting with proving a useful lamma. Lemma A.3. Let x 1 and let m be an integer. Then
(2(m - k) - 1)x* 5 nd
(1 - x)²
Proof. First, observe that
m-1
m-
(2(m- k)) 1)x* (2m - 1)
Σ
x 2 Σ
kr
k=
Since x # 1, the geometric scries partial sum is
m-1
x
1 1-
0=
17

---

Published in Transactions on Machine Learning Research (07 (025))
A standard result (or by differentiating the geometric series) gives
m I. mxx" t + (m 1)
kx
(1 - x)²
k=
Multiplying the desired inequality
(2(m - -k) 1) $ xx + 1)
(1 x²
by (1 (1-x) 2)2 (which is positive for x 1) reduces the task proving
(2m - 1)(1 - x²² (1 - x) - -2 (x mxx + (m 1) S (x 11.
A short computation shows that the left-hand side simplifies to
(2m 1) (2m + 1)x xm + m+1
so that the inequality becomes
2m 1 S (2m + 1)x.
Since x > 1, the last inequality holds, and the proof is complete.
Lemma A.4. Suppose Tm = r < T and Assumption hold, then for any m 2 0,
+ a
E Um umll [1] S
nn + O (max a.) ")),
(1+ + van)
where
21113,0² (1 - P)2 711-m(11pp) p)) pp' p)* Pn (1 + np)? (p),
ag (1 + na)2
Proof. Throughout this proof, we assume E |4 1 # 0, because if [| - - = 0 almost surely for all
j, then there is no error to bound and we directly obtain E [|um 1m = 0, which trivially satisfies the desired inequality.
We begin by introducing the following notation:
Xa.j = pt- 8Ns
P..
P
= n Σ 8
i=f+1
1=0
From Lemma A.2
the difference Um
Um can be expressed in terms of the deviations Å, A. Specifically,
faxes (A)-A) 8,,0 + a (XRs(A, - A)(8, + +8,)){
Tm
km
1=0
m
k-1
Σ
Xb. (Ag 2) A1) (X,,0 + =a,))
r7
A 0
1=
Letting
Zy m- k-1 Σ
(n A) (RA0 +ar,) c
18

---

Published in Transactions on Machine Learning Research (07/202) 2025)
we then write
m m-1
Um
mm
17
Zk-
(Bounding Blllim umll using the martingale sequence)
Zx. is the martingale sequence with respect to its
input Fk-1 (Ag...., Ak-1) because from 1B[Ax-1 A] = 0, we have
Zl F ... Fol = BlZzK F2
k-2
pm m-mk-
Σ
(A, - A) (8,0 + + + XR [A- - a] (82-100 + aP +
1=0
1k-1
Using this property, we compute
m-1
m-1
m-1 m-1
E
Σ
Σ [|]
+2 Σ
Σ
[
k=
k= j=k+1
m-1 m-1
[|| 2 +2 Σ
Σ
E
2% Biz, FN
k=
j=k+1
=2p
m-
m-11
[1| 2 +2 Σ
(m
1
K)F
1x
k=
k=
Hence,
m
E [1um
Fm]12] S 17 'E
Σ Zz
K=
m
(2(m- k)) 1)r km
(A) -A) (8,.0 + aP;) c
1)
m
k-1 (2(m- k) 2)-2(m-
|| - A) + +88,) )c
|| F
1
k=
j=0
where the last equation uses the property of martingale difference sequence:
E
[[xxs(A; A) (X,00 + aX;)) 1 Fo = O.
(Bounding the inner variance norm) By expanding the variance, for some random vector X, we have a common property
E X - || = E | 21 [xT [[[]]]] + |
=E [XX -1BBX11
So applying this to
X = = A) &V,)c,
19

---

Published in Transactions on Machine Learning Research (07 (07/2025)
condition on Fy-1 to get
[[(A - () + (8 155-11 = E|IlA, + + [ - [| + a8z)c|? 1F-11
= 11A(8
Bounding ||A,|| S from Assumption we get
E [IKA A) alj)el?] B [183,0 a
(Handling the geometric sums)
Substituting into the earlier inequality, we obtain the error bound:
k-1 Σ
E Im Fm Fmll S n ||e| 8²
Σ k=
(2(m k) 1)) -
Xxx
E [ + +
(37)
1=0
To simplify the notations for the bound in the summation, let
R R= (1 - -n) + nd, R = (1 m) + np, r= 17) + nj..
Then bound each of | ||&l and || obtained by geometric terms in R R:
|xxمl R
|x s R, | S i
P0--11
1 Nd
s=0
It follows that
E 11 at,11 s (11 (1-0
1-
?
Pe + +a-"
1
2
1 R
1 R
2
1- γ 17
R20 +2
17
1 - R
1 R
where the last inequality uses ||a + b12 < a| + 2/ exprossions:
Let's define the following symbols to simplify our
Rp
R2 = R2, Rp = R2?
Tp = 22
Tp
R
Ki
1-- 17
K2 2
7)
IT2 - R7
1 R
1 R
Rp
Using these definitions, we can continue our derivation:
2
27
X 1 |8
+ a8,1
< (2(--511)
2 α
17
R2j
17
R
ft.
R
Ki P2] K2
R/
Ki RVj. + R j-1
122
R
1
R
Ki R2
K2
R
Rh-1
KITT + K2TT
20

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Using the above derivation, 37 is expressed
m-1 (2(m- - k) I)we
1
E mm Fm 12] S n 0²
r Rp
11
22
(38)
We then split this into decaying part and non-decaying part.
First, consider the first term of the decaying factor in (38). Recalling T > 1 implied by E [1 All] # and Lemma
A.3 we have
--1
m-
R1
Rp
p
m-
Σ
(2(m k) 1)
1
Rp
m P
(2(m - k)) - 1)w (xī -1) 1)
k=
5=0
R=0
(WTT) m
1
Rp
p
(WIT) - 1)
Rg
K1
1
rp
mm11
Cp
fg
R
Rt
Rp
1
7p
Rmp
1
p
K1
r
Rip
Rp)
R
1
2
O Rp:
(39)
r4
Now consider the second term of 138
which the non-decaying part. Similarly, we use Lemma A.3 and
#2 >
from E
A|l2] # O, having
m- (2(m k:) 1)w*
1
172
(2(m -k) 1)w* (T 1)
K2
Rp
p
12
K2
Rp
I2
1=0
(WT2)) m (WTT2 + 1)
S
K2
Rp I22
(WTT2 1)² ?
1
[+ mp
1
y 1m-1
Rp
2
1 -1
1>
1 1 Tp
K2
Rp
(1 rp)
Since the denominator can be replaced with
1- R² 1 - (1-1) -==1-((1=))+ = (2 1) (1- p)) p),
(1- - ft)2 = ((1
-1 (1
= 1,²(1 +)) P)²
(1- = (1 (1-(1- (1-
(-n((1a))))) = (2 m(1 2 (1 a))2
the non-decaying part is finally expressed as
12 Rt
1
(2(m- k) 1)
Σ
T2
1 1+ - (1(1-(11-aa)))
(40)
12 - () p)) (2 -n(1 - =))
1=0
21

---

Published in Transactions on Machine Learning Research (07/2025)
By injecting the decaying term 39 and non-decaying term
440
into 38
present the final result as:
(1 - - -1) - (1-n(1nn))a
7 + (1
B [I|um Jm 1m - 20
n11 - a))2 + O (max( an) )m)
1+ 1+ (1-nn(1-a)) a))²
on + (py:
(2 11 aa)
1 0
Jn + 0 (max anl an)
(1+ + van)
A.6 Convergence Analysis of Deterministic MixedFP-KM
Lemma A.5. Suppose Tim = - and Assumption 3 hold, then for any 0,
110m d, = (,)
where
7 np)) Pn = (1
an = (1 - + ma))
Proof. From Lemma A.1
we have
- (orux owwk)
Um
7
m-
k-1
m- m-
ph p*c aapccaa
D p
7
km0
==0
h-1
7
m- k-1
ptcc σn
p
(41)
k=
s=0
Using the identity for the matrix inverse hn-1 X (I-X""
X²) x) and the fact that
- M(1-P)-1c= - (1- T/(1 - A)))
1n - d.f.
the term inside the sum in
41 can be rewritten as
k-1
pb c + an
pk-s-1 c pp on (11 - ( - P)
8=0
p (1 - A adzf.
By letting D= = 1 - a(1 -A)-1 and substituting them into
41 we obtain
m-
Fm d =
m k-1 m (p&
De ad,f) - daff
p m mm k- pt Dc
an
(1 )m)) def def
1n pk D dd
k=0
22

---

Published in Transactions on Machine Learning Research (07 (025)
Using ||a + b112 s ||a]12 + 1|b11, have
m-1
||c|12 +2 (1 pm
2
Um - dofll' 2nj |D||
Σ
p m k - Rb.
||l||
1--p
k=0
1m
pm
2
pm
21j² |||||
||c| +
||e|12
R r
1 p
O (max R, r, r)am) = (max an) mm)
A.7 Main Convergence Theorems
Theorem 4.5 (MixedFP-KM), Let a E 1), When 7m =1 < T-p and Assumplions to 3 hold, then for any m 2 0,
1 on
E [m d2S113 S
n + O (max an) )
(1 + v(n)
where
27,2.3.22
Tn = (1 () 711 - pp'
Pn= = (1 7 np)
On
an (1 17 + 1/a22
Proof. We divide the difference into
Um - d (Um - Um) + (Um dgf).
where 7m the deterministic counterpart. Then
[l d
E [| Fm Fm + EE [[|Tm d
where the cross term vanishes because E [um Um] Lemma,
0. Applying Lemma A.4 (for ||
4m II) and
A.5 for [mm. dpfll) proves the claim.
We obtain Theorem as the special case of Theorem 4.5 where rj with additional consideration for
Theorem 4.1 (MixedFP). Suppose Assumptions to hold and E [0, 1], then for any m 2 0,
01 (pT")
ifa {0,1). σ1 + O (max (p1,))p )m otherwise,
E [][mm - d2/[117] <
1+a1
+ /on
where
217
P1 === p. a1 a²
(1- 1) p2)
Proof. For a E [0, 1), this result is the special case of Theorem 4.5 where = 1, which gives an = =a² = an
213,2" Pg 8 Pl, and on [1-0) ai When 1, 6m is equivalent to êm-1 whose error is given by
Corollary B.
23

---

Published in Transactions on Machine Learning Research (07 (07/2025)
B Theoretical Analysis of StocFP with the Stochastic KM Iteration
B.1 Convergence of StocFP with the Stochastic KM Iteration
Theorem 4.6 (Stochastie KM iteration of StocFP). Suppose 1mm = < T-B and Assumplions Then for any m 2 0,
to 3 hold.
E [m fofll S on 0 (p%)
where on and Pn defined in Theorem 4.5
Proof. The proof divided into two parts: the variance bound (that is, the stochastic error between Wm and wn), and the doterministic crror bound (i.c., the bias u| - dgfl)?). We then combine these two parts.
From 29 the error in the w sequence can be written as
m-1 Σ
pm- 1m--l1(Aj - A) {
[I 1=0
II P }
Wm
1j
-A)
c
==
R=0
Nizk+1
Recall
Xms = pm-1-) 800== = II P. p -I
A.
k=0
k=t
Then,
Σ XmslAAg -1) A) (8 + +8i)e.
Wm Wm 7)
Because B[A,] A from Assumption 3 the sequence
W, == X A) (x
forms a martingale difference sequence with respect to the natural filtration F, oA,o...., Az). Hence, by orthogonality of martingale differences and the Cauchy-Schwarz incquality,
E [||m
19]
E
17
W
1=0
2)²B
Σ j=0
W
n"²
E
[||
j=0
By expanding the term in the summation, we have
E [1] = E || [ - A8 + +85)d11]
= =1IXm,jl1 2 [I( 1) + 83)cc12]
Conditioning on F5-1, we get
E [[ A)($,,,0 + +1)dd12 1 = [| + 830013 1
1|A({8, + 8,)0
24

---

Published in Transactions on Machine Learning Research (07 (07/2025)
where we used the E ||| - [[X]||]] E [||X|]]] ||[[]]|? From Assumption 3 gives us
have ||A,|| s A. which
E [I(CA A) [3)d1 s e |133,0 + 8/1²
Therefore,
E [I[W]1 [[mm 1 ||] || + Y313]
(42)
Furthermore, since the norm of P = I- nI A) is bounded by ||||| I R with R (p- p). Thus,
| || I Rm-1-j
Similarly, the product of the stochastic matrices satisfies
|| S Re,
|& S 1
R p
7) (1 RP)
1 R.
k=
with R= (1 (1-n) 7) + np 1 17(1 p). Therefore,
| + 8,113 $
R' +
(1
R¹)
1 1 R
< 1
1]
R23 +
11
1 pt
1 R
Substituting them into 42 gives
[m Wm
1 s n/ /2||c|1
0²
m2
R 2m 1-1)
2 211-1
2
A2j +2
17
(43)
1 1t.
1- R
f=0
By letting
2
2
R R
Kg =
7)
K4
1)
Rp R2,
RG nt,,
T2
R
1 R
1
R
We simplify the notation of as
43 as
E [[lwm
wm 11 S p2 12 Σ
(neTf + kATH))
k=
Since R= 1 - 1/(1 - p) < R = 17(1 - P) < from > and O <
T-p:
1
we have
k 1
k
R
(Namt + kaTh) 1---1
Ki I
+
K4
1
π2
1
k=
fm
1
K3
+ K4
Rp
Rp
1
R2p
fm
K4
Rp
+
Ap
17(1 - p2)( 17(1 p1) 0 (?
2
25

---

Published in Transactions on Machine Learning Research (07 (07/202) /025)
It follows that
E Wm Tmll2] & 22
7 (2- p)) (1 +0 ((1 (11-11) (1))
We now derive in full the bound for the deterministic error Um dgfll?. The deterministic iteration is given by
To 4
Wm+1
Püm + ne, with P= 17(1 1 A).
Unrolling this recursion yields
m
Wm pm c +
pm pm-1-k
k=
Changing the index (letting j = m 1- 1-k)
k) we have
m-
Wm = p + Σ
pi pm p P)-
5=0
Subtracting d,f from um yields
Wm dzf [pm c + 17 P
pm (1- P1- - 17(1 - P
= - 11 11(1 - pm 1] P) 7e
pm - 1/p (1 - P
pm [c-n(I PP-
Since d,f = = - P, we rewrite this as
Wm - -
f = pm (c dxf)
As |P|| < 7/11 p) =: R holds under Assumption 3 we cleduce that
m defll = O (11- 17(1 pm)
We now decompose the overall
Wm dxf (wm Wm) Ww dzf).
Since E [Wm Wm] = 0, we have
E [[lwm dz < [[]w wm + [m dz.l11
Inserting the bounds derived above yields
E [m dofl13] < 20 (1- (2 17(1 p)) (1 +0 ((1 771 2m))
We obtain the error bounds for StocEP Corollary B.1 (StocFP and StocRB). Suppose Assumptions
4 StocRB 5
as the special case of Theorem 4.6 where 1 =
to hold. Then for any m 2 0,
E [l d. d E [l]mm d, A Sa + (
where on and pr defined in Theorem 4-1
26

---

Published in Transactions on Machine Learning Research (07 (075)
B.2 Comparison with the Existing Result
Theorem 4.6 for StocPP with the stochastic KM iteration (StocFP-KM) slightly differs from bound derived from results in the original paper (Grazzi et al. 2021
Grazzi et al. 20221 Theorem 5.1) presents a non-asymptotic error bound for StocFP-KM with a decreasing stepsize, while the bound for a fixed slepsize is partially provided their Theorem 4.1, which non- asymptotic version of our Theorem [.2 By combining their Theorem 4.1 and their proof of Theorem 5.1, we can derive the following result for StocFP-KM with a constant stepsize.
Theorem B.2 (Fixed stepsize and asymptotic version of Grazzi et al. 2021 Theorem 5.1)). Let 02 2(0" +0""+ and suppose <イ S (1-0)*
1 and Assumplions to hold. Then for any m 2 0,
Itoz
Wm - d, d,fll S 2 (1 p1 p)2 + +0 ((1 - 1(11
Although both Theorem 4.6 and Theorem B.2 evaluate the same method, the derivations differ. We chose to evaluate the error of the existing method using approach similar to our own analysis in Theorem 4.5 in order to ensure fair comparison between our proposed method and that of Grazzi et al. (2021). Con- sequently, our Theorem 4.6 does not unfairly penalize the method in Grazzi et al. (2021) in fact, there are
cases where Theorem 4.6 even yields tighter bounds than Theorem B.2
One example is obtained by choosing
p p= 0.1, p 0.2, 1 10
These values satisfy the side conditions we have
0 <<p <<<<l 1, s. and n - By dividing both errors by 2nL p²
(((-M))))))
((-M))))
22 11( p))(1 +0
م))
+0
p²) )²)
[1 p²) (1 D)3
=0.844
=0.990
#1.247
This is an example where both the non-decaying term and the decaying term of our Theorem 4.6 are superior to Theorem B.2 obtained from the existing results.
C Impact of Approximate Inner Solutions on Estimation Accuracy
In the previous sections, we analyzed the convergence properties of MixodFP-KM and StocFP-KM under the assumption that the inner problem is solved exactly. However, in practical applications, we often work with approximate solutions of the inner problem. This section investigates how the accuracy of hypergradient estimation is affected when using inexact inner solutions, providing theoretical guarantees for the error bounds in such scenarios.
In short, the errors induced by inexact inner solves are consistent between the mixed and stochastic estima- tors, and the variance term, which MixedFP KM is expressly designed to curb, remains unchanged.
We suppose that Jacobians and gradients are evaluated at an approximate inner solution Rds, writing : z()), introducing smoothness of stochastic Jacobians and gradicnts.
Assumption (Smoothness with respect to the inner paramcter). There erist finite constants LA, L > such that for every Rd,,, xr, x' Rdd, and & E X,
|0, x) On (aelaa,As)))) A:E1I < LAIlE ='Il, |0,s(1, A) 0,f(x',)) S Lelkæ - x'1l.
C.1 MixedFP-KM with Inexact Inner Solution
MixedFP-KM runs with Jacobians and gradients evaluated at. ã, producing the stochastic sequence
Zm (Um, Wm, Umn) m 0, 1,...
27

---

Published in Transactions on Machine Learning Research (07 (025))
A; by its expectation) is and its deterministic counterpart (same recursion, but replacing each (, denoted žm. = (im, üm. um))
We also use the following standard resolvent identity.
Lemma C.1 (Resolvent identity). Let A, A E Rd xdz be such that I - A and I À are invertible. Then
(I- -)- (I = = A) A)(T A)-
(44)
Proof. Starting from the hight side,
(1 ( - - A)(I - A) (1 - (1) -
because Á A- (1
(1 - A). Distributing the product,
(I AI A)(I A) 1 - (I AI A) = - A)- AI-
This completes the proof.
Theorem C.2 (MixedFP-Ke) P-M with inexact inner solution). Let óz = -- -A)).. Suppose Assumptions 1 to hold and that MixedFF KM is run with Jacobians and gradients evaluated at ã, generating Un. Assume also n p-p.
1 Then for all m 2 0,
1 an
B Bllimm defl?? (1 /a,,)20 + K|00211 + ((max(pp) "m)
where Ogy Pro On are from Theorem 4.5 and
1
||c|| 1
K K=2
(1 D)²
(1 p)a
Proof. Let À = O2 A), ē= 02f(E A). We decompose the total difference:
Üm f (Im- )..)) + (um i00) (i dzf).
Here, 8m is the deterministic sequence obtained by replacing cach 0z x; (x) with Å, and i.0 -II-A)- Taking norms and expectations, we use (a b)2 < 2a² + 202 and note that Ee(ū em) O, cross terms vanish in expectation. Thus,
E [ d
< E 6m 1m + (l 800) +5x0 d+f)l2
< EE @l 8m - 1 + l1o00 - flfll + 2 8oo 11 - dxfll + +l - 0001
First, replacing by À and e by in the proof of Lemma A.4 shows
1 + an
B[l@m 1m. S
(1 a.y) 2 o + (m,,
-1 1 (I A) -ccc is handled by the resolvent identity of noting that ||-- All S P. Next, 800 dnff = -I (I- -)
Lemma C.I. which implies
-1--) = (I - A -À- A)(I - -
Hence noting ē = c) G,
(I- - 8 (I A) c= (I A)- + A) (A A)(1 A)-
A1-
Because ((I-A)- -1 - and A All |o ||
cll < La|ldz|| from Assumption 4
we get
- lē- cll
À - Alll|
Le + ||c||LA
1 +
||ØØ||.
(1 2)
1 p (1- pp)
28

---

Published in Transactions on Machine Learning Research (07 (025)
Therefore its squared norm is at most
+ ||c| 1 - d < 2 <0(==2(111) (1-p)a p)2 (1 p)* ||xx||?
which is our constant term K|
Finally, by the analysis in Lemma A.5 with the single fixed matrix A and vector é, get
| 800 ll (
Recalling that 10 dzʃll is constant,
2 0m 00 cc - d 1m0m 2011 = (maxx(Pr m
Collecting all three parts completes the proof.
C.2 StocFP-KM with Inexact Inner Solution
The Tollowing result quantifies the additional error that appears when the fixed point equation for the inner variable is solved only approximately. All assumptions stated earlier remain in force.
Let be an inexact inner solution and define
8x == ã - xAA),
Am = (,
T:= 0 A).
For fixed stepsize n (0, I-p the stochastic Krasnosel'skit Mann (KM) version of StocFP is
ip č um+1 (1 mm + (AAm um + +
0,1,...
Theorem C.3 (StocFP-KK) FMM with incxact inner solution). Let Assumptions to hold and let om be generaled above. Then for every 2 0,
Blllüm d fll S o. K||0x|1 Op""),
(45)
where on and pn
pn= = =1 + n are the constants of Theorem 4.6 and K is defined in Theorem C.
Proof- Set Ā == == (ix A) and = 0, A)- Let Wm be the deterministic KM sequence obtained by replacing every Am with A and define the limit toc (1 A) 1
Exactly as in the proof of Theorem
C.2 decompose the error and bound the squared norm:
E fl d < E @m
uml! || + [1(im wxx) + (w 2
EE [m um [ + || d 21 " do.fl + 10mm -0x0
With A replaced by À and by 7, the variance bound in Theorem 4.6 yields E| L - wm < o, + Opp)))
The resolvent identity Lemma C.1
gives Th dzf (I A) c) (I A) A)(I A Using
[I(T - 1 < (1 - -1_ A Al < Lalldul1, and ||č cll < La 18x/ yields
|| defll [ + |للاي |(2x|]
whose square bounded by K
Applying the deterministic part of Lemma A.5 to the single matrix À gives
im W ORRT).
Therefore, the cross product 2||um Wee |l|| Rxx - dgfll because the first factor decays geometrically while the sccond is constant (O())&xI))
Collecting the three contributions gives the stated bound.
29

---

Published in Transactions on Machine Learning Research (07 (275))
C.3
Discussion
The two error bounds derived above share the same struclure and, to leading order, the some magnilude. In both Theorems C.2 and C.3 the mean square error can be written as
(non-dccaying error from variance) +
(decaying crror)
blas from inexact solution o(p7) or O (maxx (P)
on (1v00,220)
where the constants K Ou, Cups and Py are identical (up to the harmless factor (1 + an)((1 + van))) in the two approaches.
The key observation is that the variance component is unaffected by inner-loop bias. As shown in the proof of Theorem C.2 the perturbation åx = ã =(A) influences the estimator only through the resolvent difference, giving rise to the deterministic term K| [ it leaves unchanged the second-moment bound of the martingale difference sequence that drives the variance. Consequently, the variance-reduction mechanism built into MixedFP KM is preserved even when the inner fixed-point equation solved only approximately.
Experiment Setup Details
In this section, we provide detailed information about the setup used in our empirical evaluations.
D.1 Effect of Mixing Rate
The experiment in Section 5.] evaluates the estimation error of hypergradients in the synthetic setting where the Jacobian matrix is a random variable. We implemented the following configuration:
The experiment used the following parameter settings: dimension of inner-parameter d, = 10; dimension of outer-parameter d, 10; number of parent matrices .I,...,4k Hr, for sampling 1000; eigenvalue bound
parameter 10-2; maximum iterations M 1000; scale parameter values y E {10-3, 3 10 1000, 2, 1_10! 1.0): and mixing rate values corresponds to StocRB.
E (0, 0,000001,.01,0000 0.01, 0.1, where a 0 corresponds to StocFP and 1.0
Any matrix in Hi...., ".....,H, Ha was constructed to ensure that each eigenvalue follows a uniform distribution in [0, e]. The matrix H; E (HI,..., Hx} was constructed through the following process: first, generating a random matrix R Rd, with entrics sampled from standard normal distribution; then creating a symmetric matrix S = (R+ (R+RT)/2 )/2; next, computing the QR decomposition of S to obtain orthogonal eigen- vectors Q; followed by sampling eigenvalues { uniformly from [0, 1 e]; and finally constructing
H. Q diag(A)...., QT.
The coelficients E Rd=, d E Rdd, and B E Rdaxds were sampled independently from a uniform distribution in [0, 1] for each element.
D.2 Compare with Existing Approaches
In this section, we provide detailed information about the setup and results in Section 5.2
D.2.1 Setups
The synthetic problem setting differs from that in Section 5.1 in the following aspects: we used fixed y 1.0 instead of varying it over [0, 1], the dimensions of inner-parameter and outer-parameter À were both set to 100, we ran experiments with 10 different random seeds (instcad of 100), and we used 10,000 iterations for each method. For the hyperparameter optimization task, we used binary classification with the Adult Income dataset. For the influence estimation task, we standardized the Fashion MNIST data with zero mean and unit variance. For the meta learning task, we used the min-max scaling on the California Housing dataset with the regularization parameter 8 0.1. Table 2 shows the dataset sizes and the values of y for f(x, A) (T +F)x x BA used in the three weal-world tasks.
30

---

Published in Transactions on Machine Learning Research (07/2025)
Task
Dataset
Ptrain
Rval
dg
dx
7
Hyperparameter Optimization
Adult Income
5000
5000
14
14
100
Influence Estimation
Fashion MNIST
5000
5000
784
5000
10-2
Meta Learning
California Housing
5000
5000
8
8
10-1
Table 2: Experiment settings for the real-world tasks.
D.2.2 Influence Estimation
We define the influence of the j-th instance A as the change in validation loss [xx A) after excluding the j-th instance from the training dataset with samples, denoted by
4yf f(x(11)))))))_)11) 14)) - f(x(1),1)
where 14y {0, 1j is vector with ones except for denote an index of the instance in the dataset and g(x, x; A; {i) with (5d)ig represent the inner loss of the i-th instance & masked by the i-th element of x. Then, the inner solutions x(1) and x(1y;) are defined as the minimizers of an with X = 1 and 14j, respectively:
at the j-th index. Let i Uniform([1,.... NN)
x(1) arg arg min E, [g(2, 1;
x(14)) arg min E, [g(2, 1z;;
o(1+;) can be seen as the solution of the inner problems in which the j-th instance is excluded from the dataset, since Eggg(z, 14y; i)] ト Ei= g(2, 1) 2/4 E 6, ((;;)() where dG is the Kronecker delta. Under this formulation, the j-th element of the hypergradient - -d2, (x((), at À == provides a linear approximation of the influence of the j-th instance, because
f(x(14)), 144) - f(x(1), dxf( (1 +)+ 0 0[1 -113)
= = =-d,,(=(W)))))= 0(1).
A similar definition of the influence and its approximation were also presented in Koh & Liang (2017) Recalling æ(x, A) N E, Ong(2, ); at X 1, we have
O.p = H, 82 -r,(=))))(1rE)))
where H H= = E, (G(rE) From
2a
we have
- - T -
which is the downweight version of their Eq. (2).
For the initialization of manee-paramtteraa x, used samples from the normal distributions scaled by constant factor for both logistic regression models used in the hyperparameter optimization and influence estimation tasks and linear regression models for the meta learning. For the outer-parameter X, the initialization varies by task. In the synthetic task, is initialized with values sampled uniformly from [0,1]. In the Hyperparameter optimization task, serves as a regularization coefficient initialized with a small constant value. In the influence ostimation task, represents loss weights for each training sample, which is a ones- vector. In the meta learning task, À functions as biased regularization term to the initial value of the model's parameters
Any inner-problem optimization was performed using the Adam optimizer with learning rate of 0.01. To rule out the effect incurred by inexact x(A), for any task, we used the full-batch inner loss to compute gradients for Adam and ran 1,000 epochs to ensure the convergence.
We optimized the hyperparameters using grid search find the configuration that minimizes the aver- age error over the last M/10 iterations. For each method, we performed grid search over the following hyperparameter ranges:
31

---

Published in Transactions on Machine Learning Research (07/2025)
10-3, 3.1 10 10-2, Mixing rate 1 - 10,10-4,2x 10 100$,10 10
10-2, x2100_
Stepsize 7, E (10-3,2 x 10 5, 10
10 2, 2 10 10 3,5x
10 .2x 10 .5 10 1.0}
For decreasing stepsize schedule: 3 E 20, (10, 20, 50, 100, 200, 500, 2000, 5000, 10000} with the constraint 3 م 8
20, 50, 100, 200, 500, 1000, 2000, 5000, 10000) and ó e
D.3 Results
For completeness, we provide additional experimental results that complement the main findings presented in Section 5.2
Fig. g presents the same experimental results as Fig. 2 but with error bars indicating the standard deviation
across multiple runs. The error bars represent the standard deviation above and below the mean squared error across runs with different sequences (Em)
We also provide computational cost comparison since our MixedFP and MixedFP-KM methods require computing two vector-Jacobian products (VJPs) per iteration,
Fig. | shows the results from Fig. and Fig. with the horizontal axis changed from algorithm iterations to the number of VJP computations. This transformation stretches the horizontal axis for MixedFP and MixedFP-KM by factor of two, since our approach requires computing two VJPs per iteration: A, üm and um as shown in 66b.. This represents worst-case analysis assuming VJP computation dominates the total cost and that two VJPs require exactly twice the computation time.
The results show that our methods' advantage is reduced when evaluated iteration count, though performance remains competitive with existing approaches. Note that this represents a worst-case analysis assuming VJP computation dominates total cost and our two VJPs require exactly twice the computation time. In practice, the actual overhead is less than 2x because the two VJP computations share commnn computations (particularly the forward pass, as they use the same samples). Moreover, the two VJP computations can be parallelized, potentially hiding the wall-clock impact.
wall-clock basis compared to
32

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Synthetic Problem
Hyperparam. Optimization
10a
10-2
f2 error from V,f
10
10-2
109
10
2k
4k
6k
8k
10k
0k
2k
4k
6k
8k
10k
Influence Estimation
Meta Learning
12 error from V,f
Copttd ha & fiv aen
10-2
10-4
10-3
0k
10k
2k
4k
6k
8k
10k
Number of iterations
Number of iterations
StocFP
Kt (StocFP KM-Iteration with constant n)
MixedF++KM- (MixedFP-KM with constant n)
StocRB
L (StocFP KM-iteration with decreasing
MixedFP+KM (MixedFP-KM with decreasing n)
MixedFP
Figure 3: Results corresponding Fig. 1 with added visualization of variance across runs. Each line shows the mean squared error over 10 trials with different. Jacobian sample sequences (Am)men, with error bars indicating standard deviation.
33

---

Published in Transactions on Machine Learning Research (07 (07/2025)
Synthetic Problem
Hyperparam. Optimization
10a
10-1
Va16
12 error from V,f
101
10-2
10a
10-3
4k
6k
8k
10k
0k
2k
4k
6k
8k
10k
Influence Estimation
Meta Learning
f2 error from Vaf
10-2
10-4
10-3
0k
10k
2k
6k
10k
Number of VJP computations
Number VJP computations
StocFP
StocFP+KM- (StocFP KM-iteration with constant n)
MixedFP (MixedFP-KM with constant
StocRB
StocFP+KM FStocFP KM-iteration with decreasing
MixedFP+KMI (MixedFP-KM with decreasing n)
MixedFP
Figure 4: Comparison stochastic hypergradient estimation methods evaluated by number vector- Jacobian products (VJPs). The horizontal axis represents the number of VJPs, effectively stretching MixedFP and MixedFP-KM curves by factor of two. This provides worst-case analysis of wall-clock performance, assuming VJP computation dominates total cost and our two VJPs require exactly twice the computation time. Method-specific parameters were re-optimized under this new metric.
34

---

Published in Transactions on Machine Learning Research (07 (07/2025)
E
Additional Experiments: Hyperparameter Optimization
In this section, we evaluate the practical impact of improved hypergradient estimation accuracy actual bilevel optimization performance. While Section 5.2 demonstrated that our MixedFP-KM method achieves superior hypergradient estimation accuracy, it remains unclear whether this cheoretical advantage translates to better optimization performance in real bilevel optimization tasks. To address this question, we conduct bilevel optimization experiment using the hyperparameter optimization setting adopted in Section 5.2
E.1 Settings
We employ the same inner and outer problem formulation as the hyperparameter optimization task in Section 5.2 where the outer-parameter represents regularization coefficients for each dimension of the The inner problem minimizes the binary cross-entropy loss on the training dataset, while the outer problem minimizes the validation loss with respect to the regularization paramcters. In addition to the training and validation splits used in Section 5.2 we introduce separate test set of 5,000 samples to evaluate the final model performance after the outer optimization.
mener-parameter
We configure the bilevel optimization with 100 outer oplimization steps using SGD with a learning rate of 20.0, and 100 inner optimization steps per outer iteration using Adam with a learning rate of 0.01. The hypergradient estimation depth set to 100 iterations (reduced from 10,000 Section 5.2 for computational efficiency), and we conduct experiments across five random seeds to get the averate test accuracy. We compare all hypergradient estimation methods from Section 5.2 StocFP, StocRB, MixedFP, StocFP+KM StocFP+XM4, MixedFP+KM+, and MixedFP+KM4, with all method-specifie hyperparameters (n, a, B, and d) set to the optimal values determined by grid search as in Section 5.2
E.2 Results and Discussion
Fig. presents the learning curves for both validation and test losses across all methods. Methods incorporat- 'sann-Maan ing the stochastic Irainone iteration (MixedFP+KML, MixedFP+KM+, and StocFP+KMJ) demonstrate substantially better performance than their non-KM counterparts (StocRB, StocFP, and MixedFP). Among these stochastic-KM-bassd methods, MiXedFP+KMI exhibits the steepest initial loss reduction and the highest test accuracy, which correlated with the hypergradient estimation accuracy results shown in Fig. G where MixedFP+KMI achieved the smallest error. These results validate that the improved hypergradient estimation accuracy demonstrated in Section 5.2 can translate practical benefits in bilevel optimization.
35

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
x.4 10
0.80
5.2 x10
5 10
Test accuracy
0.79
Test loss
4.8
0.78
4,6 10
0.77
10-2
0.76
4.2 x 10-1
0.75
20
60
100
0
60
80
100
Number of outer steps
Number outer steps
StocFP
KK (StocFP KM-Iteration with constant n)
MixedFP+KM-1 (MixedFP-KM with constant MixedFP+KM (MixedFP-KM with decreasing n)
StocRB
KtocFP+KM (StocFP KM-iteration with decreasing
MixedFP
Figure 5: Test loss and accuracy curves for hyperparameter optimization of regularization coefficients across different hypergradient estimation methods. Curves show the average values over 5 runs with standard deviation error bars.
f2 error from V2
10-2
StocFP
StocRB
MixedFP
StocFP+KM-> (StocFP KM-iteration with constant n)
StocFP+KM1 L (StocFP KM-iteration with decreasing n)
MixedFP+KM- (MixedFP-KM) with constant
MixedFPKKMI (MixedFP-KM with decreasing n)
10
Number of iterations
Figure 6: Hypergradient estimation errors at the initial outer-step using the optimal method-specific param- eters (n,, a, 3, and 8) determined by grid search, which are then used in the bilevel optimization experiment shown in Fig.
36

---

Published in Transactions on Machine Learning Research (07 (07/)025)
def mixed_fpp ____p__mm_____ km
f: Callable [ [Tensor, Tensor], Tensor],
phi: [allable[[Tensor, Tensor, Tensor], Tensor],
Tensor,
lambda_: Tensor,
samples: Sequence [Tensor],
etas: Sequence
alpha: float, M: int
Tensor:
Compute outer gradients and initialize v, w, u torch.zeros. _eros_like($ (x)
(c,) = torch. autograd grad (x, lambda_), x)
w, = c.clone(), c. c.
for m in range (M)
eta etas [m] samp samples [m]
phi wi with fixed lambda
phi fixed = lambda x_: phi (x_, lambda_, samp)
Compute VJPs
(Au,) torch. autograd functional. "jp(phi_fixel (x,), (H,))
autograd.fuucttonall
(Au,) aurchauttgrad functional. happppi_fieed (x,), (u,))
Update v, w,
v (1 eta) * eta * (alpha * (v w + (1 alpha) w)
w (1 eta) * w eta * (Aw + c)
u (1 eta) * eta
* (Au)
return detetcch()
Figure 7: Minimal Python implementation of the MixedFP-KM, corresponding 9
F Python Implementation of MixedFP-KM
The following Python pscudocode Fig. 7 provides explicit implementation of the MixedFP-Ki algorithm, as introduced and discussed in Section 4.2
In this implementation, the vector-Jacobian product for the stochastic mapping Åm (x), is computed using the automatic differentiation tools PyTorch. The itcrative updates for the variables (v, w. u) and other notations correspond to those given
37


Reviews:
[{"review_id": "bYm8b9tpSQ", "reviewer": "Reviewer_StT7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- The paper is extremely well written. From the choice of the notations to the comparison side by side of the algorithms, the authors managed to provide a very clear manuscript to an idea that is actually quite technical.\n- The application may be slightly niche, but the idea is original and the results demonstrate clearly the advantage of their method on several benchmarks.\n- The theory gives interesting insights that are well illustrated with carefully crafted experiments. \n\n**Weaknesses**\n- A snippet of code could be beneficial for the paper to grab adoption. As the authors mention, the implementation requires typically hvps (or vjps of the implicit equation) that readers could be unfamiliar with.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a new estimation technique of hypergradients from stochastic estimates that achieves variance reduction over previous techniques by mixing them. The authors detail clearly two previous methods and how they propose to combine them. In addition, they propose to blend their method in Krasnosel’ski˘ı-Mann iterations to further enhance it. They provide theoretical bounds on the expected mean squared error showing improvements over previous methods. They provide numerical experiments showcasing the advantages of the method that their theory forecasted.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Question:\n- \"The error comprises two distinct decay rates, and for sufficiently large m, the slower decay rate eventually dominates.\" I don't understand: for a given problem $\\rho_1$ and $\\alpha$ are fixed so there is no reason for the decay rate to switch from one decay to the other.\n- \"While both StocFP (4) and StocRB (5) provide the same hypergradient approximation in expectation, they yield different biases\". Biases are generally the difference between the expected estimator and the true estimate. So if they are the same in expectation, I would expect them to have same biases. Did you want to speak about e.g. their variance? Sorry if I'm missing something obvious.\n\nTypos:\n- $\\sigma_1$ is used both Theorem 4.1 and Theorem 4.2, consider using two different notations (same for $\\sigma_2$\n- Maybe use $\\eta_m$ instead of $\\eta_k$ in the sums of Theorem 4.2\n- In Lemma 4.3, say \"the function $F(z)=...$\", avoid using math symbols as entities of a phrase.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "JdH5M0G1BX", "reviewer": "Reviewer_8hMc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**:\n- The idea to construct a mixed estimator is natural and carefully executed.\n- The mathematical derivations are sound and based on well-established tools from stochastic approximation and fixed-point theory.\n- The theoretical analysis precisely quantifies the error bounds and clearly identifies when the proposed method offers improvements.\n\n**Weaknesses**:\n- The improvement in variance bounds is modest.\n- While MixedFP shares the same order of computational complexity as StocFP and StocRB, it requires updating and maintaining multiple iterates at each step, resulting in increased per-iteration computational cost.\n- The theoretical analysis for decreasing stepsizes is limited to asymptotic convergence, without non-asymptotic error bounds, while decreasing stepsizes are used in some experiments.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes MixedFP, a stochastic hypergradient estimation method that combines StocFP and StocRB updates to reduce the variance in hypergradient estimation. The authors further apply the stochastic Krasnosel’skiĭ-Mann iteration to MixedFP to guarantee almost sure convergence to the true hypergradient. Theoretical analyses provide variance bounds, and empirical results on synthetic and real-world tasks support the approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Clearly acknowledge that MixedFP requires updating multiple iterates per step, resulting in increased per-iteration computational cost compared to StocFP and StocRB.\n- Explicitly state that the theoretical analysis for decreasing stepsizes is limited to asymptotic convergence, without non-asymptotic error bounds, and reflect this limitation earlier in the main text.\n\nAlso, some suggestions :\n- Report experimental comparisons based on wall-clock time or normalized computational effort, rather than iteration count alone, to ensure a fair evaluation.\n- Discuss that the experimental setting assumes full-batch convergence of the inner problem, which does not reflect the more challenging stochastic or approximate scenarios where hypergradient estimation is typically critical.\n- Briefly comment on the practical difficulty of satisfying the unbiased Jacobian estimation assumption in real-world large-scale applications.", "broader_impact_concerns": "Not applicable.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "hmuYlBBb0Z", "reviewer": "Reviewer_M2Gy", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n* The proposed hypergradient estimator is a simple modification of known estimator which achieves a more accurate and faster hypergradient approximation. This improvement is performed without any supplementary oracle computation.\n\n* Numerical experiments corroborate the theoretical findings.\n\n* The estimator comes with sound theoretical guarantees.\n\n### Weaknesses\n* **Clarity and precision**: Clarity has to be improved, in particular by clever notation choices (see requested changes for details) and better wording.\n\n* **Approximation in $x(\\lambda)$**: In the hypergradient approximation problem, there are two bottlenecks: the resolution of the linear system (which is treated in the paper) and the approximation of $x(\\lambda)$. However, the approximation of $x(\\lambda)$ influences the hypergradient precision, also in the linear system since the exact hypergradient uses $(I-\\partial_x \\varphi(\\textcolor{red}{x(\\lambda)}, \\lambda))^{-1}\\partial_x f(\\textcolor{red}{x(\\lambda)} \\lambda)$ while the approximate hypergradient uses an approximation of $(I-\\partial_x \\varphi(\\textcolor{red}{\\tilde x}, \\lambda))^{-1}\\partial_x f(\\textcolor{red}{\\tilde x}, \\lambda)$ where $\\tilde x$ is an approximation of $x(\\lambda)$. This approximation is not accounted for in the analysis while being fundamental in hypergradient estimation.\n\n* **One-shot estimation**: Related to the previous weakness, the paper only studies one-shot estimation of the hypergradient. For instance, for the different experiments, comparing the performances of the whole pipeline when this hypergradient estimator is plugged would have been interesting (that is, not only the hypergradient error but also for example the performance in optimization for the hyperparameter optimization task). \n\n* **Reproducibility**: The source code of the experiment is not provided, degrading the reproducibility.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper considers the problem of hypergradient estimation. It proposes MixedFP, a combination of the known algorithms StocFP and StocRB, which provably improves upon both algorithms with a faster convergence. Then, the authors show that the MixedFP algorithm can is actually the iterations of a contracting fixed point map. This enables them to improve MixedFP by proposing MixedFP-KM, which implements Krasnosel’skii-Mann iterations. Non-asymptotic $l^2$ error bounds are provided.\nNumerical experiments of hypergradient estimation are provided on a synthetic problem, the hyperparameter optimization problem, the influence estimation problem, and meta-learning.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "### Clarity/precision\n* **Confusion between outer cost $f$ and value function $\\lambda\\mapsto f(x(\\lambda), \\lambda)$**: The clarity could be improved by a clever choice of notation. Indeed, denoting $\\nabla_\\lambda f$ the hypergradient is confusing since it can be confused with the gradient of $f$ with respect to $\\lambda$. Giving a name to the value function, that is, the function $\\lambda \\mapsto f(x(\\lambda), \\lambda)$ as classically done in the literature, would help to clarify the exposition. Moreover, page 3, the sentence \"The hypergradient is the gradient of an outer-cost $f(x(\\lambda), \\lambda)\\in\\mathbb{R}$ with respect to the outer-parameter $\\lambda\\in\\mathbb{R}^{d_\\lambda}$, where the inner-parameter $x(\\lambda)\\in\\mathbb{R}^{d_x}$\" is really confusing, since the hypergradient is NOT the gradient of the outer function evaluated at $(x(\\lambda), \\lambda)$ but it is the gradient of *the composition of the outer function and the solution map of the inner fixed point problem*.\n\n* **Implicit differentiation and approximate implicit differentiation**: Sometimes the paper calls \"Approximate implicit differentiation\" things that correspond to \"Implicit differentiation\" (without \"approximate\"). Implicit differentiation provides a closed-form formula of the hypergradient using the implicit function theorem. Approximate implicit differentiation is a set of techniques that compute an approximation of the exact hypergradient by replacing in the hypergradient formula, the exact value of the inner variable $x(\\lambda)$ and the exact value of the linear system solution $(I-\\partial_x \\varphi(x(\\lambda), \\lambda))^{-1}\\partial_x f(x(\\lambda), \\lambda)$ by approximation. For this reason, I think that the sentences \"A method called approximate implicit differentiation (Pedregosa, 2016; Lorraine et al., 2019; Rajeswaran et al., 2019) considers the following form derived from the chain rule and implicit function theorem...\" and \"Approximate implicit differentiation (Lorraine et al., 2019; Pedregosa, 2016) uses the implicit function theorem and the chain rule to rewrite the hypergradient in a form that involves an inverse matrix...\" should be rephrased.\n\n* **Hessian-vector product**: Hessian-vector  (HVPs) are inappropriately mentioned on pages 4 and 6 for the following reasons:\n    1. The paper considers the general setting where the inner problem is a fixed point iteration and not necessarily an optimization problem. As a consequence, the hypergradient computation in the different algorithms does not necessarily involve HVPs but involves Jacobian-vector products.\n    2. Page 4:  \"$\\hat A_mw_m\\in\\mathbb{R}^{d_x}$ can be calculated in $O(d_x)$ time using the Hessian-vector product technique\". The term \"Hessian-vector product\" only refers to the product between an Hessian matrix and a vector, without presuming the way it is computed. They can be computed naively by first computing the Hessian matrix and then multiplying it by the vector, or more efficiently using automatic differentiation techniques [1, 2]. Moreover, the $O(d_x)$ should be justified by a reference.\n\nThus, the sentences evoking HVPs should be rephrased in light of these two points.\n\n* **Influence estimation experiment**: I don't get the setting of the influence estimation experiment. In particular, I don't see why the hypergradient catches the behavior of the performance when a sample is excluded. Moreover, I don't find the formulation of the problem proposed by the paper in the provided references [4, 5]. I think this part needs to be rewritten.\n \n* **Page 3**: \"which is equivalent to solving (2b) in Section 3.1 by stochastic gradient descent\". It should be precised that, since $\\varphi$ is contractive, solving (2b) is equivalent to minimizing a strongly convex quadratic function.\n\n* **Page 3**: \"$x(\\lambda)$ is defined by the stationary point of a mapping...\" -> \"$x(\\lambda)$ is defined as the fixed point of a mapping...\"\n* **Page 5**: \"Section 4.2 explains that unbiased hypergradient estimation is achieved\". I find this sentence misleading, since, due to a finite number of iterations, the $v_m$ is not an unbiased estimator of $\\nabla_x f$. However, it is an asymptotically unbiased estimator. \n\n### Experiments\n* **Error bars**: the experiment results should report the variability of the results over the different trials.\n\n* **Numerical comparison of the algorithms**: the complexity of an iteration of each algorithm is not identical. Therefore, comparing them with respect to wall-clock time rather than iteration number would be fairer.\n\n* **Experiments with the full pipeline**: as mentioned in the weakness section, an experiment, for instance with hypergradient descent for hyperparameter selection, would enhance the paper.\n\n\n### Appendix\nI found some typos and some steps that should be a bit more detailed in the appendix, without affecting the final result.\n\n* **Page 13, proof of lemma 4.3**: It seems that there is a confusion between $F$ and $T$.\n* **Theorem 4.4**: If the $\\sigma_2$ of the paper is the same as $\\sigma_2$ in [3], then a square is missing on the denominator.\n* **Appendix A3.**: Precise explicitly that the fixed step-size case is considered in this section to avoid confusion.\n* **Lemma A.1**: Maybe I missed it but I did not find the definition of $\\bar{\\alpha}$ (which I presume is equal to $1-\\alpha$).\n* **Induction in Lemma A.1**: At the induction step, (19) is assumed for a given $m\\geq 0$. Then it is shown that $v_{m+1} = \\eta \\sum_{k=0}^m r^{m-k}(\\alpha u_k + \\bar\\alpha w_k)$ which is nothing else than (19) which is assumed to be true. To be rigorous, the induction step should show that, given $(19)$, we have $v_{m+2} = \\eta \\sum_{k=0}^{m+1} r^{m+1-k}(\\alpha u_k + \\bar\\alpha w_k)$. This remark holds also for (20) and (21).\n* **Page 18**: Why $\\lVert\\hat A_j - A\\rVert^2\\leq \\hat\\rho^2$?\n* **Equation (38)**: Should $\\lVert \\hat A_j - A\\rVert^2$ be removed in the right-hand-side?\n* **Page 19**: It is not obvious to me why the equality holds $\\mathbb{E}\\left[\\lVert \\hat X_{j,0} + \\bar\\alpha \\hat Y_j\\rVert^2\\right] = \\left(\\left(1-\\bar\\alpha \\frac\\eta{1-\\hat R}\\right)\\hat R^j + \\bar\\alpha \\frac\\eta{1-\\hat R}\\right)^2$, thus it deserves a bit more details.\n* **Page 19**: $\\bar \\alpha \\frac\\eta{1-R} = (1-\\alpha) \\frac1{1-\\hat \\rho}$: $R$ should be replaced by $\\hat R$.\n* **Pages 20 and 21**: Contrary to what is written, the inverse matrix identity is useless to get $\\eta(I-P)^{-1}c = \\nabla_x f$ since it is straightforward from the definition of $P$. But it is useful after when handling the term $P^kc + \\bar\\alpha\\eta \\sum_{s=0}^{k-1} P^{k-s-1}c$. Thus the sentence introducing the identity should be moved to the correct place.\n* **Page 21, first equation**: $c$ is missing in the r.h.s. of the first line.\n\n\n[1] Barak A Pearlmutter. Fast exact multiplication by the hessian. In *Neural computation*. 1994.\n\n[2] A. Griewank and A. Walther. *Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation*, Second Edition. Society for Industrial and Applied Mathematics, 2008.\n\n[3] R. Grazzi, M. Pontil, and S. Salzo. *Convergence properties of stochastic hypergradients*. In AISTATS, 2021\n\n[4] P. Koh and P. Liang. *Understanding black-box predictions via influence functions*. In IICML, 2017.\n\n[5] R. Khanna, B. Kim, J. Ghosh, and S. Koyejo. *Interpreting black box predictions using fisher kernels*. In AISTATS, pp. 2019.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)