{"id":"335738d1-4f59-4a6c-a804-baf7e0115c45","arxiv_id":"2411.18100","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":7.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"A zeroth-order proximal-gradient method with Gaussian smoothing is proved to find approximate stationary points of stochastic bilevel inverse problems with non-smooth convex lower levels, at oracle complexity O(epsilon^-3) in the nonconvex case and with explicit accuracy in the convex case.","lead":"This paper develops a derivative-free stochastic optimization algorithm for bilevel inverse problems, where one learns regularization or measurement-design parameters without computing gradients of the reconstruction map. It gives formal complexity guarantees and demonstrates the approach on signal denoising and angle selection for X-ray tomography.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The inexact lower-level complexity rests on Lemma 4.10, whose Appendix B.4 proof assumes independence between ξ1 and ξ2 and a global Lipschitz F that the quadratic loss in Sec. 8 does not satisfy; Assumptions 5/6 are also unverified for that example.","rationale":"The reader identified Assumptions 5(S.2)/6 as the weakest load-bearing assumption. I agree those assumptions are unverified and can fail in the small-regularization regime, and the paper's bounded-box numerics do not rescue the statements on all of Y. However, my stress-test found a more immediate internal obstruction: the inexact lower-level analysis, which is a central advertised contribution and the only regime implemented numerically, depends on Lemma 4.10. Its proof factorizes expectations of dependent random variables and requires a global Lipschitz constant for the squared-error loss, which is not available on R^{nx}. This is not merely a missing verification but a proof step that is invalid for the paper's own Section 8 loss. Because the flaw is repairable by adding explicit boundedness/independence or by replacing the global Lipschitz assumption with a localized/moment-based condition, the appropriate verdict remains CONDITIONAL rather than REJECT: the exact-case complexity seems coherent, but the paper must be revised to state and prove the inexact results under assumptions that cover its numerical experiments.","tokens_in":40569,"tokens_out":10771,"duration_ms":105958,"concrete_test":"Re-derive Lemma 4.10 for F(x,ξ1)=||x-X||^2 and D=KX+Z without factoring conditional expectations. Using Hölder instead of the invalid product factorization, show that E||a_{k+1}||_*^2 is bounded by a constant times β_k^2/η^2 only if quantities such as E[||xβ-x*||^{2p/(p-2)}] or E[||X||^{4p/(p-2)}] are finite. Then check these moments for the Section 8.1 model: if D has Gaussian noise and X has bounded support, the fourth-moment condition may hold, but the global-Lipschitz assumption on F still fails; if no such moment bound can be established from Assumptions 3-6 alone, Lemma 4.10 is false as stated and the inexact complexity bounds require an additional assumption.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The most load-bearing concern is not the exact-case O(ε^{-3}) bound but the inexact lower-level analysis. Lemma 4.10 is the only step controlling the bias terms a_{k+1}, b_{k+1} in decomposition (4.19); without it, Theorems 5.6 and 6.1 and the inexact-number complexity statements have no proof. Appendix B.4 derives E[||a_{k+1}||_*^2 | F_k] ≤ C_F β_k^2/η^2 by factoring E[lip0(F)^2 ||xβ-x*||^2 ||U||^2] into a product of marginal expectations. This factorization implicitly requires ξ1 and ξ2 to be independent conditionally on F_k, but in the paper's own data model D = KX + Z, ξ1 = X and ξ2 = D are dependent. The same proof also assumes a global Lipschitz constant lip0(F(·,ξ1)) with finite L2 norm, while the implemented loss F(x,ξ1) = ||x - X||^2 in Sec. 8.1 is not globally Lipschitz on X = R^{nx}; no boundedness of X or truncation of the loss is stated. Thus the inexact oracle model, which is exactly what the numerical sections use, is not covered by the theorems as written. The reader's Assumptions 5(S.2)/6 concern is real, but it is secondary: even if the lower-level solution map is Lipschitz, the bias-variance chain of Lemma 4.10 fails for the quadratic loss on an unbounded domain without additional assumptions.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops a zeroth-order (derivative-free) stochastic proximal-gradient method for bilevel optimization problems arising in inverse problems, where the lower-level problem is a convex, possibly nonsmooth variational problem and the upper-level objective is an expected risk over regularization or design parameters. The method uses Gaussian smoothing and finite-difference estimators of the hypergradient, and treats both exact and inexact lower-level solutions. The paper proves iteration and oracle complexity bounds for the prox-gradient mapping in the nonconvex case, an optimality-gap bound in the convex case, a Goldstein-stationarity complexity statement, and reports numerical experiments on one-dimensional signal denoising and on Radon-based experimental design.","tokens_in":40935,"tokens_out":7097,"duration_ms":64548,"significance":"If correct, this is a useful contribution to derivative-free stochastic bilevel optimization under nonsmooth lower-level problems. The paper gives explicit convergence rates, handles inexact lower-level solves as a bias term, and connects the smoothed complexity analysis to a Goldstein-type stationarity criterion. The exact-case analysis is largely standard and plausibly correct, and the numerical experiments illustrate the intended application. However, several load-bearing assumptions are not satisfied by the paper's own numerical model, and the inexact-bias analysis contains an unjustified independence factorization. These issues affect the theoretical claims that cover the reported experiments.","major_comments":[{"comment":"The proof of Lemma 4.10 factors conditional expectations such as E[lip0(F(·,ξ1))^2 ||xβ−x*||^2 ||U||^2 | F_k] into a product of marginal expectations after invoking \"independence of the triple (ξ1, ξ2, U)\". The paper does not assume ξ1 and ξ2 are independent; indeed Example 1.1 and Section 8 take ξ1=X and ξ2=D=KX+Z, which are dependent. The same factorization appears in the proof of Lemma 4.9. Since Lemma 4.10 is the only control of the bias terms a_{k+1}, b_{k+1} in decomposition (4.19), Theorems 5.6 and 6.1 are not proved for the paper's own data model. This is load-bearing because Section 8 implements exactly the inexact estimator (4.9).","section":"§4.7, Lemma 4.10 and Appendix B.4"},{"comment":"Assumption 3 requires F(·, ξ1) to be globally Lipschitz with finite first moment of the Lipschitz modulus. The upper-level loss used in both experiments is F(x, X)=||x−X||^2 on X=R^{n_x}, which is not globally Lipschitz. No truncation or boundedness of X is stated in Section 8. Consequently the variance bound in Lemma 4.7, the estimate (4.3), and all subsequent complexity results do not apply to the numerical setting as written. The authors should either state a modified assumption (e.g., bounded X with a truncated loss) or prove the needed bounds under the quadratic loss plus boundedness of the iterates.","section":"§8.1 and Assumption 3"},{"comment":"Assumption 5(S.2) requires y↦x*(y, ξ2) to be Lipschitz uniformly in ξ2, and Assumption 6 requires square-integrability of the Lipschitz modulus of H(·, ξ). The paper does not establish either property for the lower-level problem (8.1) with TV and Tikhonov regularization, and these properties are not consequences of the stated µ-strong convexity. In particular, as λ or τ approach zero the solution map can become highly sensitive, and the theorems are stated on all of Y while the experiments constrain y to a bounded box. The hypotheses of the convergence theory are therefore not verified for the numerical example.","section":"Assumptions 5 and 6; §8.1"},{"comment":"When bounding the b_{k+1} term, the proof uses E[||U||^{2p/(p−2)}]^{(p−2)/p} ≤ n. For a standard Gaussian vector this is false for p>2; Lemma A.1 yields an upper bound of order n + p/(p−2) after raising to the appropriate power. The resulting constant C_F should depend on p and n. This is repairable and does not change the rates, but as printed the lemma is not correct.","section":"Appendix B.4, bound following eq. (4.20)"}],"minor_comments":[{"comment":"The proof contains a repeated typo: the displayed bound for ||tilde G_{η,α_k}(y_k) − G_{η,α_k}(y_k)|| ends with \"= ||∆W_{k+1}||_*^2\", but it should be ||∆W_{k+1}||_*.","section":"Appendix B.2, proof of Lemma 5.1"},{"comment":"The sentence \"the regularization parameters λ, ν and the smoothing parameter ν\" repeats ν; it should presumably read λ, τ and ν.","section":"§8.1.1"},{"comment":"The notation is inconsistent: the lower-level regularizer is r2(x) in (1.3), but two paragraphs below the text refers to \"a regularizer/penalty function r2(y)\" in the upper level; this should be r1(y).","section":"§3, after eq. (1.3)"},{"comment":"The numerical sections state that Algorithm 1 is run with inexact lower-level solutions, but Definition 4.6 requires an L^p error bound over the sample. The paper should explain how the implemented stopping criterion enforces the L^p condition, or state that the experiments are heuristic and not covered by the inexact theory.","section":"§8.2"}],"recommendation":"major_revision","confidential_remarks":"The exact-case analysis is a reasonable contribution and the paper fits the journal's scope. The main revision priority is the inexact-oracle section: repair the independence factorization in Lemma 4.10, correct the Gaussian moment bound, and either verify or explicitly restrict the assumptions for the quadratic-loss experiments. I would be willing to re-review a revised version."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Worth a serious referee: the paper is a genuine extension of zeroth-order stochastic MPEC analysis to a proximal composite setting with inexact lower-level solves, and the Goldstein-subdifferential transfer is clean. The exact-case complexity results are new and appear correct. The paper deserves a major-revision decision, not desk rejection.\n\nThe soft spot is the inexact lower-level analysis. Lemma 4.10 is the only control on the bias terms a_{k+1}, b_{k+1} in the error decomposition (4.19). Its proof in Appendix B.4 factors conditional expectations of lip0(F(·,ξ1))^2 and ||x^β−x^*||^2 into a product, which is only valid if ξ1 and ξ2 are separately independent conditional on the filtration. The paper never states this, and its own data model (Example 1.1 and Section 8.1) sets ξ1=X and ξ2=D=KX+Z, so they are dependent. That breaks the proof.\n\nSecond, the same proof and Assumption 3 require F(·,ξ1) to be globally Lipschitz with an integrable modulus. The quadratic loss ||x−X||^2 used in the experiments is not globally Lipschitz on R^n. The parameter y is restricted to a box in the numerics, but the theorems are stated on all of Y, and no boundedness or truncation of X is added. So the inexact oracle model that the numerical section actually uses is not covered by the theorems as written.\n\nThere are smaller, fixable issues: the bound on E||U||^{2p/(p−2)} in Lemma 4.10 is numerically off (n should be roughly sqrt(n+2p/(p−2))), and the oracle count in Corollary 5.4 appears to miss a factor of two per function evaluation. These do not change the structure.\n\nWhat holds up: the exact-case analysis is careful, the smoothing-to-Goldstein-stationarity argument is a real contribution, and the numerical evidence shows the method works in practice on the tested problems. The assumptions 5(S.2) and 6 are strong and unverified for the examples, but that is secondary to the Lemma 4.10 gap.\n\nWho is it for: researchers in stochastic bilevel optimization and derivative-free methods for inverse problems. I would send it to a serious referee, with the expectation that the inexact analysis be repaired or clearly separated from the experiments. My own verdict would be major revision.","headline":"Original exact-case theory, but the inexact lower-level bias analysis relies on an unstated independence of ξ1 and ξ2 that the paper's own data model violates.","tokens_in":41492,"tokens_out":5690,"would_cite":false,"duration_ms":47852,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["90C15","90C26","90C56","65K10"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that regularization parameters and measurement designs in inverse problems can be learned by a stochastic bilevel method that only queries objective values, with total oracle and lower-level complexity O(ε^{-3}) to reach…","keywords":["stochastic bilevel optimization","zeroth-order methods","derivative-free optimization","inverse problems","regularization parameter learning","Gaussian smoothing","Goldstein subdifferential","complexity analysis"],"falsifier":"Compute the Lipschitz modulus of the solution map $x^*(y,\\xi_2)$ for the lower-level problem (8.1) as the Tikhonov parameter $\\lambda\\to 0$ at fixed data; the theory requires $\\|\\mathrm{lip}_0(H(\\cdot,\\xi))\\|_2<\\infty$. If data and parameters exist with $\\lambda\\to 0$ where this modulus is unbounded or not square-integrable, then the variance bound $s^2=(4+n)^2\\|\\mathrm{lip}_0(H(\\cdot,\\xi))\\|_2^2/m$ is infinite and the $O(\\varepsilon^{-3})$ guarantee collapses; a direct numerical check measures the empirical variance of $V_{\\eta,m}(y)$ along such a sequence and observes it diverging as $\\lambda\\to 0$.","tokens_in":40330,"feed_emoji":"🎯","tokens_out":9915,"duration_ms":78501,"temperature":0.7,"pith_summary":"The paper proposes a stochastic bilevel optimization method for learning regularization parameters (and, in optimal experimental design, measurement policies) in inverse problems from data, when the lower-level problem is a convex nonsmooth variational reconstruction. Its central claim is that this task can be solved without gradients: the algorithm estimates the gradient of a Gaussian-smoothed version of the hyperobjective using only function values, allows the lower-level problem to be solved only inexactly, and still provably converges. The proof supplies finite-time complexity bounds: $O(\\varepsilon^{-2})$ iterations and $O(\\varepsilon^{-3})$ total calls to the function oracle and lower-level solver to reach small expected prox-gradient norm in the nonconvex case, an explicit convergence rate for the convex case, and a relaxed Goldstein stationarity guarantee for the nonconvex case. If correct, this gives a disciplined, oracle-cheap alternative to hyperparameter grids and manual tuning in data-driven inverse problems.","feed_headline":"Gradient-free bilevel learning in O(ε^{-3}) oracle calls","feed_subtitle":"A function-value-only algorithm with inexact inner solves provably learns regularization parameters for inverse problems.","key_machinery":"The load-bearing object is the Gaussian-smoothed hyperobjective $h_\\eta(y)=\\mathbb{E}[h(y+\\eta U)]$ with $U$ standard Gaussian, whose gradient $\\nabla h_\\eta(y)=\\mathbb{E}[(H(y+\\eta U,\\xi)-H(y,\\xi))/\\eta\\, BU]$ can be sampled without derivatives. Finite differencing in random Gaussian directions gives the zeroth-order estimator $V_{\\eta,m}$, whose variance is controlled by the square-integrable random Lipschitz modulus of $H(\\cdot,\\xi)$ (Lemma 4.7). The prox-gradient mapping $G_{\\eta,t}(y)=(y-\\mathrm{prox}_{t r_1}(y-t\\nabla h_\\eta(y)))/t$ serves as the stationarity measure, and its stochastic analogue is what Algorithm 1 drives to zero. Proposition 4.5 converts smallness of $\\nabla h_\\eta$ into membership in a Goldstein $\\delta$-subdifferential of the original nonsmooth $h$, which is what turns the complexity bounds into a meaningful relaxed-stationarity guarantee.","core_discovery":"The central discovery is a derivative-free proximal-gradient scheme (Algorithm 1) for the hyperobjective $\\Psi(y)=\\mathbb{E}[F(x^*(y,\\xi_2),\\xi_1)]+r_1(y)$, where $x^*(y,\\xi_2)$ is the unique solution of the lower-level variational problem. The method replaces the unavailable hypergradient by a multi-point finite-difference estimator built from Gaussian random directions and i.i.d. data samples; Gaussian smoothing turns the nonsmooth expected hyperobjective into a $C^{1,1}$ function $h_\\eta$, so the prox-gradient machinery applies. With step sizes $\\alpha_k=2\\beta/(\\mathrm{lip}_1(h_\\eta)\\sqrt{k})$ and batch sizes $m_{k+1}=a\\sqrt{k}$, the paper proves $\\mathbb{E}\\|G_{\\eta,\\alpha_\\kappa}(y_\\kappa)\\|^2\\le\\varepsilon$ after $O(\\varepsilon^{-2})$ iterations and $O(\\varepsilon^{-3})$ function and lower-level oracle calls (Corollary 5.4), and the same rate with controlled inexact lower-level solves (Corollary 5.7). In the convex case it bounds the expected optimality gap explicitly in terms of the smoothing parameter $\\eta$ and the inexactness schedule (Theorem 6.1), and in the nonconvex case it shows the output is an $(\\varepsilon_1,\\varepsilon_2)$-stationary point in the Goldstein sense with $O(\\varepsilon_1^{-4}\\varepsilon_2^{-2})$ iteration complexity (Theorem 7.1).","pith_inferences":["A consequence the authors leave implicit: the theory requires the Lipschitz modulus of the lower-level solution map to be square-integrable over the whole parameter space, which is not verified for the TV/Tikhonov lower-level problem; restricting $y$ to a compact box bounded away from zero, as the experiments do, is a practical but unproven workaround.","Because the variance constant grows with the dimension $n$ of the parameter space $Y$, the method is most attractive when few hyperparameters are learned, such as a handful of regularization parameters or a softmax-parametrized measurement policy.","A testable extension is to replace i.i.d. Gaussian directions with quasi-Monte Carlo or multilevel Monte Carlo sampling, which the paper flags as future work; the complexity bounds here predict the logarithmic factor would shrink if the cumulative batch sum were reduced.","The complexity transfer to Goldstein stationarity suggests the same algorithm could handle other nonsmooth composite upper-level problems, such as group-sparse hyperparameter learning, provided the same smoothing and square-integrability conditions hold."],"forward_implications":["If the bounds hold, users can learn regularization parameters and experimental designs with no hypergradient computation: only objective-function samples and an inner solver for the lower-level problem are needed.","For the recommended decaying step-size and growing batch-size schedule, reaching expected squared prox-gradient norm $\\varepsilon$ requires $O(\\varepsilon^{-2})$ iterations and $O(\\varepsilon^{-3})$ total inner solves.","In the convex case the expected optimality gap vanishes at a rate explicit in $\\eta$ and the $\\beta_k$ schedule, letting the user trade inner-solve accuracy against total cost.","In the nonconvex case the output is certified as a relaxed Goldstein stationary point, a necessary relaxation because no first-order method can certify ordinary $\\varepsilon$-stationarity in finite time.","Constant step sizes and constant batch sizes converge only to a plateau, so the schedule is essential for the guarantee."],"supporting_citations":[{"why":"Supplies the Gaussian smoothing lemmas that make the estimator unbiased, bound the smoothing error, and give the delta-subdifferential inclusion used throughout.","marker":"[40]"},{"why":"Provides the key transfer result (Proposition 4.5) connecting smoothed gradients to Goldstein subdifferentials, used for relaxed stationarity.","marker":"[35]"},{"why":"Introduces the stochastic zeroth-order method for mathematical programs under equilibrium constraints that this paper extends to nonsmooth proximal composite bilevel problems.","marker":"[13]"},{"why":"Supplies the inexact derivative-free bilevel learning setting and the lower-level strong convexity constants and numerical setup used in Section 8.","marker":"[18]"},{"why":"Defines the Goldstein subdifferential used for the relaxed stationarity notion.","marker":"[26]"},{"why":"Shows no first-order method can find ordinary epsilon-stationary points of nonsmooth nonconvex functions in finite time, motivating the relaxed Goldstein criterion.","marker":"[50]"},{"why":"Analyzes inexact hypergradients in bilevel learning and underpins the inexact lower-level oracle model.","marker":"[19]"}],"fun_headline_variants":["Gradient-free bilevel solver hits O(ε^-3) oracle calls","Function-value-only bilevel learning at O(ε^-3)","No-grad bilevel optimization for inverse problems","Bilevel inverse problem learning, gradient-free, O(ε^-3)"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The proof hinges on the lower-level solution changing Lipschitz continuously with the learned parameters, with the random rate of change having finite second moment; this can break down when the regularization parameter approaches zero, which is precisely the regime of greatest interest in inverse problems.","fun_headline_variants_meta":{"raw":{"variants":["Gradient-free bilevel solver hits O(ε^-3) oracle calls","Function-value-only bilevel learning at O(ε^-3)","No-grad bilevel optimization for inverse problems","Bilevel inverse problem learning, gradient-free, O(ε^-3)"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00068,"raw_usage":{"total_tokens":3151,"prompt_tokens":1069,"completion_tokens":2082,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":685,"completion_tokens_details":{"reasoning_tokens":2019}},"tokens_in":685,"tokens_out":2082,"duration_ms":15633,"temperature":1.0,"reasoning_tokens":2019,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T11:30:58.926437+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the Lipschitz modulus of the solution map $x^*(y,\\xi_2)$ for the lower-level problem (8.1) as the Tikhonov parameter $\\lambda\\to 0$ at fixed data; the theory requires $\\|\\mathrm{lip}_0(H(\\cdot,\\xi))\\|_2<\\infty$. If data and parameters exist with $\\lambda\\to 0$ where this modulus is unbounded or not square-integrable, then the variance bound $s^2=(4+n)^2\\|\\mathrm{lip}_0(H(\\cdot,\\xi))\\|_2^2/m$ is infinite and the $O(\\varepsilon^{-3})$ guarantee collapses; a direct numerical check measures the empirical variance of $V_{\\eta,m}(y)$ along such a sequence and observes it diverging as $\\lambda\\to 0$.","supporting_citations":[{"cited_title":"Random gradient -free minimization of convex functions","cited_arxiv_id":null,"evidence_quote":"Supplies the Gaussian smoothing lemmas that make the estimator unbiased, bound the smoothing error, and give the delta-subdifferential inclusion used throughout."},{"cited_title":"Shanbhag, and Farzad Yousefian","cited_arxiv_id":null,"evidence_quote":"Introduces the stochastic zeroth-order method for mathematical programs under equilibrium constraints that this paper extends to nonsmooth proximal composite bilevel problems."},{"cited_title":"Ehrhardt and Lindon Roberts","cited_arxiv_id":null,"evidence_quote":"Supplies the inexact derivative-free bilevel learning setting and the lower-level strong convexity constants and numerical setup used in Section 8."},{"cited_title":"Optimization of lipschitz continuo us functions","cited_arxiv_id":null,"evidence_quote":"Defines the Goldstein subdifferential used for the relaxed stationarity notion."},{"cited_title":"Complexity of finding stationary points of nonconvex nonsmooth functions","cited_arxiv_id":null,"evidence_quote":"Shows no first-order method can find ordinary epsilon-stationary points of nonsmooth nonconvex functions in finite time, motivating the relaxed Goldstein criterion."},{"cited_title":"Analyzing inex act hypergradients for bilevel learning","cited_arxiv_id":null,"evidence_quote":"Analyzes inexact hypergradients in bilevel learning and underpins the inexact lower-level oracle model."}],"review_version":1}