Paper deep dive

Adaptive Robust Estimator for Multi-Agent Reinforcement Learning

Zhongyi Li, Wan Tian, Jingyu Chen, Kangyao Huang, Huiming Zhang, Hui Yang, Tao Ren, Jinyang Jiang, Yijie Peng, Yikun Ban, Fuzhen Zhuang

Year: 2026Venue: arXiv preprintArea: cs.AIType: PreprintEmbeddings: 121

Intelligence

Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 96%

Last extracted: 3/26/2026, 2:30:31 AM

Summary

The paper introduces a robust multi-agent reinforcement learning (MARL) framework for collaborative reasoning, addressing interaction-level ambiguity and reward noise. It proposes two main components: Dual-Agent Answer-Critique-Rewrite (DACR), a structured interaction protocol for credit assignment, and Adaptive Robust Estimator (ARE), a method for robust batch experience mean estimation that stabilizes policy optimization under heavy-tailed and noisy reward signals.

Entities (5)

ARE · algorithm · 99%DACR · framework-component · 99%MARL · field-of-study · 98%GRPO · baseline-algorithm · 95%Qwen2-VL-2B · model · 92%

Relation Signals (3)

DACR → ispartof → Robust MARL Framework

confidence 98% · consisting of two components: Dual-Agent Answer-Critique-Rewrite (DACR) and an Adaptive Robust Estimator (ARE).

ARE → improves → Policy Optimization

confidence 97% · ARE provides robust estimation of batch experience means during multi-agent policy optimization.

DACR → enables → Credit Assignment

confidence 95% · enabling explicit attribution of each agent's marginal contribution to its partner's performance.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Multi-agent collaboration has emerged as a powerful paradigm for enhancing the reasoning capabilities of large language models, yet it suffers from interaction-level ambiguity that blurs generation, critique, and revision, making credit assignment across agents difficult. Moreover, policy optimization in this setting is vulnerable to heavy-tailed and noisy rewards, which can bias advantage estimation and trigger unstable or even divergent training. To address both issues, we propose a robust multi-agent reinforcement learning framework for collaborative reasoning, consisting of two components: Dual-Agent Answer-Critique-Rewrite (DACR) and an Adaptive Robust Estimator (ARE). DACR decomposes reasoning into a structured three-stage pipeline: answer, critique, and rewrite, while enabling explicit attribution of each agent's marginal contribution to its partner's performance. ARE provides robust estimation of batch experience means during multi-agent policy optimization. Across mathematical reasoning and embodied intelligence benchmarks, even under noisy rewards, our method consistently outperforms the baseline in both homogeneous and heterogeneous settings. These results indicate stronger robustness to reward noise and more stable training dynamics, effectively preventing optimization failures caused by noisy reward signals.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

120,214 characters extracted from source content.

Expand or collapse full text

Adaptive Robust Estimator for Multi-Agent Reinforcement Learning Zhongyi Li* Wan Tian* Jingyu Chen Kangyao Huang Huiming Zhang Hui Yang Tao Ren Jinyang Jiang Yijie Peng† Yikun Ban† Fuzhen Zhuang† Abstract Multi-agent collaboration has emerged as a powerful paradigm for enhancing the reasoning capabilities of large language models, yet it suffers from interaction-level ambiguity that blurs generation, critique, and revision, making credit assignment across agents difficult. Moreover, policy optimization in this setting is vulnerable to heavy-tailed and noisy rewards, which can bias advantage estimation and trigger unstable or even divergent training. To address both issues, we propose a robust multi-agent reinforcement learning framework for collaborative reasoning, consisting of two components: Dual-Agent Answer–Critique–Rewrite (DACR) and an Adaptive Robust Estimator (ARE). DACR decomposes reasoning into a structured three-stage pipeline—answer, critique, and rewrite—while enabling explicit attribution of each agent’s marginal contribution to its partner’s performance. ARE provides robust estimation of batch experience means during multi-agent policy optimization. Across mathematical reasoning and embodied intelligence benchmarks, even under noisy rewards, our method consistently outperforms the baseline in both homogeneous and heterogeneous settings. These results indicate stronger robustness to reward noise and more stable training dynamics, effectively preventing optimization failures caused by noisy reward signals. Our code is available at https://github.com/bhai114/ARE. Machine Learning, ICML 1 Introduction Figure 1: The scoring results of the reward model. We randomly selected questions and answers from the MATH500 and Gaokao2023 datasets, and then used the large model to generate the answers. The first and third figures show the scoring results of the reward model for the answers generated by the large model. The second and fourth figures are the corresponding moment generating functions. It can be seen that the data distribution of the scoring by the reward model is biased. Large language models (LLMs) have achieved impressive performance on mathematical and logical reasoning tasks (Zhang et al., 2026; Yang et al., 2026; Liu et al., 2025a; Wang et al., 2025; Liu et al., 2025c; Huang et al., 2026a). In this line of work, reinforcement learning (RL) has emerged as an increasingly common approach, providing a compelling alternative to supervised imitation learning (Ren et al., 2025; Shao et al., 2025; Zou et al., 2025). In parallel, multi-agent collaboration has emerged as a practical way to mitigate single-agent limitations such as confirmation bias and weak self-critique: pairing complementary roles (e.g., solver vs. evaluator) or aggregating multiple attempts can surface hidden errors and improve reliability (Wan et al., 2025; Lin et al., 2025; Liao et al., 2025; Chen et al., 2025c). These trends naturally motivate training collaborative LLM systems via multi-agent reinforcement learning (MARL). However, in practice, policy optimization for multi-agent reasoning is often brittle, and stability issues can dominate the training dynamics. A major source of brittleness is the reward signal. Reasoning rewards are frequently produced by verifiers, reward models, or heuristic graders, and can be noisy, biased, and heavy-tailed, with occasional extreme outliers (Figure 1). Standard GRPO (Group Relative Policy Optimization) (Shao et al., 2024) normalizes rewards using the batch empirical mean and variance. Under heavy-tailed or contaminated rewards, a small number of anomalous samples can significantly skew batch statistics, distort advantage scaling, and trigger oscillatory updates or even training divergence (Yang et al., 2026). This issue becomes more pronounced in multi-agent settings, where interaction can introduce additional variability and non-stationarity into the return distribution (Wang & Xu, 2025; Huang et al., 2026b). To address these two challenges, we propose a robust MARL framework that stabilizes policy optimization under noisy, heavy-tailed rewards. The framework consists of two components. First, DACR is an interaction protocol that decomposes collaboration into explicit stages: each agent produces an initial answer, critiques its partner’s answer, and then rewrites its own solution conditioned on the critique. This structure decouples generation from evaluation and yields richer training trajectories than unstructured dialogues. We further adopt a lightweight cross-stage reward design to encourage effective interaction, without requiring complex credit assignment. Second, ARE is an adaptive robust estimator for MARL policy optimization. It replaces fragile batch-mean normalization of experience returns with a robust, adaptive location estimator, thereby improving training stability in the presence of outliers. We evaluate the effectiveness of our framework on mathematical reasoning benchmarks, using both homogeneous and heterogeneous agent pairs across model families and scales. We further demonstrate the generality of ARE beyond language reasoning by applying it to an aerial vision-and-language navigation (VLN) task. Across both tasks, our approach consistently improves performance over GRPO-based baselines. Overall, our contributions are: • We introduce DACR, a structured collaboration protocol with a cross-improvement reward that attributes each agent’s contribution to its collaborator’s reasoning quality, enabling more interpretable and learnable credit assignment. • We propose ARE, which replaces GRPO’s fragile batch-mean advantage normalization with a robust estimator, mitigating instability in advantage estimation under noisy and heavy-tailed rewards. • We systematically evaluate both homogeneous and heterogeneous multi-agent systems on large-scale mathematical reasoning benchmarks, demonstrating consistent improvements over single-agent and multi-agent baselines in accuracy, stability, and generalization. • We further validate ARE in an embodied aerial VLN setting with a Qwen2-VL-2B VLM, achieving higher training rewards and improved NE/SR/OSR/SPL across all splits, with the largest gains on unseen environments. 2 Methodology In this section, we present our robust MARL framework, which consists of two components: DACR and ARE. DACR centers on a structured interaction protocol that facilitates collaborative error detection, while ARE provides a robust optimization procedure that mitigates the heavy-tailed rewards inherent to multi-agent settings. Figure 2: The proposed DACR interaction framework. Two agents, π1 _1 and π2 _2, act as both solvers and peer-evaluators. The reward mechanism attributes credit based on standalone accuracy and the marginal improvement induced in the partner’s revised solution. 2.1 Dual-Agent Reasoning Interaction Protocol In this section, we introduce DACR, whose workflow is illustrated in the left panel of Figure 2. DACR decouples the reasoning process into a dual-agent system in which two policies, π1θ _1^θ and π2ϕ _2^φ, alternate between generative and evaluative roles. Given a problem q, the interaction proceeds through a structured three-stage trajectory: 1. Independent Answering: Both agents generate initial reasoning chains and answers: a1∼π1(q)a_1 _1(q) and a2∼π2(q)a_2 _2(q). 2. Mutual Critique: Agents swap their answers. Agent 1 generates a critique identifying potential errors in a2a_2 (c1∼π1(q,a2)c_1 _1(q,a_2)), while Agent 2 critiques a1a_1 (c2∼π2(q,a1)c_2 _2(q,a_1)). 3. Revision: Each agent produces a refined answer rir_i by conditioning on the peer feedback received: r1∼π1(q,a1,d2)r_1 _1(q,a_1,d_2) and r2∼π2(q,a2,c1)r_2 _2(q,a_2,c_1). This protocol ensures that evaluation signals are external to the original reasoning trace, providing a more objective and adversarial diagnostic environment for solving complex problems. 2.2 Cross-Evaluation Reward Mechanism In this section, we describe the multi-agent reward design for DACR, as illustrated in the middle panel of Figure 2. To optimize DACR, we move beyond binary correctness and adopt a lightweight reward mechanism that supports fine-grained credit attribution. Let S(y,q)S(y,q) denote a scoring function that assesses the logical and numerical validity of a response y to problem q. For each agent i∈1,2i∈\1,2\, we assign stage-wise rewards for the answering and rewriting phases: Rans,i=S(ai,q),Rrw,i=S(ri,q).R_ans,i=S(a_i,q), R_rw,i=S(r_i,q). Crucially, we introduce a cross-improvement term Δi _i that credits an agent when its critique leads to a measurable improvement in its partner’s rewritten answer: Δ1=Rrw,2−Rans,2,Δ2=Rrw,1−Rans,1. _1=R_rw,2-R_ans,2, _2=R_rw,1-R_ans,1. The total return for agent i is then Ritotal=Rrw,i+γΔi,R_i^total=R_rw,i+γ\, _i, where γ>0γ>0 trades off individual solution quality against collaborative utility. This objective encourages agents to provide actionable critiques while maintaining high standards for their own final answers. 2.3 Robust Multi-Agent Policy Optimization In this section, we instantiate ARE (introduced in Section 3) for robust multi-agent policy optimization. The key idea is to stabilize advantage estimation by replacing the fragile batch mean with our proposed ARE-based location estimator. Given a batch of returns Rkk=1B\R_k\_k=1^B, we compute a robust location estimate μ~B μ_B using ARE. The resulting robust advantage is defined as AkARE=Rktotal−μ~BσB+ϵ.A_k^ARE= R_k^total- μ_B _B+ε. Each agent is then optimized by maximizing the clipped surrogate objective ℒ(θ) (θ) =[min⁡(ηAARE,clip(η,1−ϵ,1+ϵ)AARE)] =E [ (η\,A^ARE,\ clip(η,1-ε,1+ε)\,A^ARE ) ] −βKL(πθ∥πref), -β\,D_KL( _θ\,\|\, _ref), where η denotes the importance sampling ratio. By anchoring advantage estimation at a robust center, robust policy optimization suppresses gradient spikes induced by outlier rewards, yielding more stable convergence in challenging multi-agent reasoning landscapes. 3 Adaptive Robust Estimator In this section, we introduce our proposed ARE method, which can be viewed as a principled refinement of the classical Median-of-Means (MoM) estimator (Sun, 2021). MoM partitions samples into blocks, computes a mean within each block, and aggregates the blockwise estimates via a median, providing robustness to heavy tails and a minority of corrupted blocks (see Appendix B for details). However, MoM enforces robustness only at the aggregation stage: each blockwise estimate is still a sample mean and thus remains sensitive to within-block outliers. ARE addresses this weakness by replacing the blockwise mean with a robust estimator obtained by minimizing an adaptive loss (Barron, 2019), thereby enforcing robustness at both intra-block and inter-block levels and improving stability under heavy-tailed contamination. We defer the adaptive loss definition and its optimization properties to Appendix D. Specifically, following the MoM template, we partition the samples X1,…,XnX_1,…,X_n into k disjoint blocks B1,…,BkB_1,…,B_k. For each block, we compute a robust location estimate by minimizing the adaptive loss: X~j∈arg⁡minx∈ℝ⁡1|Bj|∑i∈Bjρ(Xi−x;α,c),j=1,…,k. X_j∈ _x 1|B_j| _i∈ B_jρ\! (X_i-x;α,c ), j=1,…,k. (1) We then define ARE as the median of these blockwise estimates: μ~=median⁡(X~1,…,X~k). μ=median( X_1,…, X_k). (2) The main computational challenge is solving (1), since ρ(⋅;α,c)ρ(·;α,c) can be nonconvex and involves robustness parameters (α,c)(α,c). Without loss of generality, consider a single block of size m: X~∈arg⁡minx∈ℝ⁡1m∑i=1mρ(Xi−x;α,c). X∈ _x 1m _i=1^mρ\! (X_i-x;α,c ). (3) Tuning (α,c)(α,c) via cross-validation or Lepski’s method is computationally prohibitive, especially in multi-agent settings. We therefore use alternating optimization: initialize X~(0) X^(0) and (α(0),c(0))(α^(0),c^(0)) (e.g., X~(0) X^(0) as the sample median, α(0)=1α^(0)=1, c(0)=1c^(0)=1), and iterate α(t+1),c(t+1) α^(t+1),c^(t+1) =arg⁡minα∈(−∞,2],c>0⁡1m∑i=1mρ(Xi−X~(t);α,c), = _α∈(-∞,2],\,c>0 1m _i=1^mρ\! (X_i- X^(t);α,c ), (4) X~(t+1) X^(t+1) =arg⁡minx∈ℝ⁡1m∑i=1mρ(Xi−x;α(t+1),c(t+1)), = _x 1m _i=1^mρ\! (X_i-x;α^(t+1),c^(t+1) ), (5) until convergence. Section 4 details efficient solvers for (4) and (5). 4 Optimization Perspective In this section, we outline tailored solution procedures for the nonconvex subproblems (4) and (5), respectively. 4.1 Estimation of Robust and Scale Parameters Prior to estimating the adaptive loss parameters, we first derive the following key insight by differentiating the adaptive loss with respect to α. Lemma 4.1. Fix ϵ∈ℝε and c>0c>0. For α≠0,2α≠ 0,2. Then ∂αρ(ϵ;α,c)≥ 0, ∂αρ(ε;α,c)\ ≥\ 0, and the inequality is strict whenever X≠0X≠ 0. Figure 3: The probability distribution and its negative log-likelihood corresponding to the adaptive loss function. The left panel shows the negative log-likelihood function −log⁡p(ϵ;α,1)=ρ(ϵ;α,1)+log⁡Z(α)- p(ε;α,1)=ρ(ε;α,1)+ Z(α) for different shape parameters α=0,0.5,1,2,4,∞α=0,0.5,1,2,4,∞; The right panel shows the corresponding probability density function p(ϵ;α,1)=exp⁡−ρ(ϵ;α,1)/Z(α)p(ε;α,1)= \-ρ(ε;α,1)\/Z(α). Lemma 4.1 indicates that, although it may seem natural to estimate the parameters by directly solving (4), this approach is in fact ill-posed: the objective admits a degenerate solution that pushes the shape parameter α to excessively small values. To avoid this collapse, we follow Barron (2019); Jung et al. (2024) and adopt a likelihood-based calibration. Concretely, we consider the exponential-family density p(ϵ;α,c)=exp⁡(−ρ(ϵ;α,c))/cZ(α)p(ε;α,c)= (-ρ(ε;α,c))/cZ(α) with Z(α)=∫−∞exp⁡(−ρ(u;α,1))uZ(α)= _-∞^∞ (-ρ(u;α,1))du so that Z(α)Z(α) normalizes the distribution. The corresponding negative log-likelihood is −log⁡p(ϵ;α,c)=ρ(ϵ;α,c)+log⁡c+log⁡Z(α).- p(ε;α,c)\;=\;ρ(ε;α,c)+ c+ Z(α). Figure 3 shows the induced density and its negative log-likelihood, clarifying the effect of α. Decreasing α flattens the tails, reducing the penalty on large residuals (outliers) while increasing the relative penalty near the origin (inliers). Estimation therefore involves a trade-off: smaller α discounts large errors but penalizes small ones more, whereas larger α approaches least squares, favoring inliers while charging more for outliers. This trade-off encourages the method to adapt its robustness to the empirical residual distribution, avoiding the degenerate behavior implied by Lemma 4.1. Following the aforementioned discussion, we estimate the parameters α and c by minimizing the negative log-likelihood function based on the sample, which is defined as α(1),c(1) α^(1),\,c^(1) =arg⁡minα,c>0∑i=1mρ(Xi−X~(0);α,c) = _α,\,c>0 _i=1^mρ (X_i- X^(0);α,c ) (6) +m(log⁡c+log⁡Z(α)). +m ( c+ Z(α) ). It is worth noting that in optimization problem (6), we need to jointly optimize α and c, which differs from methods such as Jung et al. (2024); Chebrolu et al. (2021) that optimize only α, and may entail greater computational challenges. 4.2 Estimation of the Mean Upon obtaining α(1)α^(1) and c(1)c^(1), directly solving (5) remains nontrivial since the adaptive robust loss is generally nonconvex. We therefore adopt the graduated nonconvexity (GNC) paradigm (Black & Rangarajan, 1996) within an iteratively reweighted least-squares (IRLS) framework (Chebrolu et al., 2021). Let ϵi=Xi−X _i=X_i-X. By matching the gradients of (5) and the weighted least-squares objective, (5) is equivalent to X~(1)=arg⁡minX∈ℝ⁡12m∑i=1mwi(Xi−X)2, X^(1)\;=\; _X 12m _i=1^mw_i\,(X_i-X)^2, (7) where the weights are induced by the loss through wi=1ϵi∂ρ(ϵi;α(1),c(1))∂ϵi.w_i\;=\; 1 _i ∂ρ( _i;α^(1),c^(1))∂ _i. However, if wiw_i are treated as unconstrained free variables and optimized jointly with X in (7), the objective admits a degenerate minimum at wi=0w_i=0 for all i (independent of X), rendering the solution meaningless. To prevent this degeneration, we introduce an outlier-process regularizer on the weights via the Black–Rangarajan duality (Black & Rangarajan, 1996). Define ϕ(z)≔ρ(cz;α,c)φ(z) ρ(c z;α,c). For the adaptive robust loss family, ϕ′(z)φ (z) satisfies limz→0ϕ′(z)=12,limz→∞ϕ′(z)=0,ϕ′(z)<0, _z→ 0φ (z)= 12, _z→∞φ (z)=0,φ (z)<0, which guarantees the existence of an analytical outlier-process function Φρ(⋅) _ρ(·) such that ρ(ϵ;α,c)=minw∈(0,1]⁡12wϵ2+Φρ(w).ρ(ε;α,c)\;=\; _w∈(0,1] 12wε^2+ _ρ(w). Consequently, (7) can be rewritten as the following constrained and regularized problem: X~(1)=arg⁡minX∈ℝ,wi∈(0,1]⁡12m∑i=1m(wi(Xi−X)2+Φρ(wi)). X^(1)\;=\; _X ,\,w_i∈(0,1] 12m _i=1^m (w_i(X_i-X)^2+ _ρ(w_i) ). (8) When initialization is poor, directly optimizing (8) remains challenging due to nonconvexity. We therefore solve it by GNC: we start from a convex quadratic surrogate and gradually introduce nonconvexity until the original adaptive loss is recovered. Concretely, we employ a shape mapping f(β,α(1))f(β,α^(1)) that interpolates from f=2f=2 (quadratic) to f=α(1)f=α^(1) (target adaptive regime), controlled by a continuation parameter β. Details are provided in Appendix E. 5 Theoretical Properties We organize this section as follows. We first study the mean estimator with a fixed shape parameter α∈(0,1]α∈(0,1] and a sample-size–dependent scale cmc_m. We then establish two complementary regimes: (i) a finite-variance regime, in which the estimator is consistent and asymptotically normal at the classical m m rate; and (i) a heavy-tailed regime, where only a (1+ϵ)(1+ε)-moment is assumed and the estimator attains a high-probability deviation bound of order (log⁡(1/δ)/m)ϵ/(1+ϵ) ( (1/δ)/m )^ε/(1+ε) for the single-block estimator in (3). Finally, we study the asymptotic theory of the median-of-means aggregation—namely, the ARE in (2)—based on the collection of block means. Throughout, the same estimating equation underlies all results; the difference between regimes is driven by how cmc_m is scaled relative to the stochastic fluctuations of the empirical score. Let X1,…,XmX_1,…,X_m be i.i.d. real-valued random variables. The goal is to estimate the mean μ:=[X1]μ:=E[X_1] (assumed to exist). Fix a shape parameter α∈(0,1]α∈(0,1] and let κα:=2−α∈[1,2) _α:=2-α∈[1,2). We use the adaptive robust loss family ρ(ε;α,c):=κα[(1+(ε/c)2κα)α/2−1],c>0,ρ( ;α,c):= _α [ (1+ ( /c)^2 _α )^α/2-1 ], c>0, with score (influence) function ψ(ε;α,c):=∂ερ(ε;α,c)=εc2(1+(ε/c)2κα)α/2−1.ψ( ;α,c):= ∂ ρ( ;α,c)= c^2 (1+ ( /c)^2 _α )^α/2-1. Heuristically, ρ is locally quadratic while ψ downweights large residuals when |ε|| | exceeds the scale c. Given a scale cm>0c_m>0 (possibly depending on m), define the empirical risk R^m(x):=1m∑i=1mρ(Xi−x;α,cm), R_m(x):= 1m _i=1^mρ(X_i-x;α,c_m), and the estimating equation g^m(x):=dxR^m(x)=−1m∑i=1mψ(Xi−x;α,cm). g_m(x):= ddx R_m(x)=- 1m _i=1^mψ(X_i-x;α,c_m). Let Rm(x):=[ρ(X−x;α,cm)]R_m(x):=E[ρ(X-x;α,c_m)] and gm(x):=Rm′(x)=−[ψ(X−x;α,cm)]g_m(x):=R_m (x)=-E[ψ(X-x;α,c_m)]. Our proofs compare g^m g_m to its population counterpart gmg_m on a neighborhood of μ to localize a stable root. For α<1α<1, g^m(x)=0 g_m(x)=0 may have multiple roots due to non-convexity. We therefore define a central root in a controlled interval. Let μ^pil μ_pil be any pilot estimator and choose a radius rm>0r_m>0 such that rm→∞r_m→∞ and rm=o(cm)r_m=o(c_m). Define ℐm:=[μ^pil−rm,μ^pil+rm].I_m:=[ μ_pil-r_m,\ μ_pil+r_m]. Definition 5.1 (Central-root estimator). If the equation g^m(x)=0 g_m(x)=0 has a unique solution in ℐmI_m, denote it by X~m X_m. The only role of ℐmI_m is to rule out spurious distant roots while keeping μ inside with high probability. 5.1 Finite-Variance Regime Assumption 5.2 (Finite variance). [X1]=μE[X_1]=μ and [(X1−μ)2]=σ2∈(0,∞)E[(X_1-μ)^2]=σ^2∈(0,∞). We use a purely asymptotic deterministic tuning: for a fixed γ>0γ>0, let cm:=m1/2+γ,rm:=mγ/4.c_m:=m^1/2+γ, r_m:=m^γ/4. With this choice, cmc_m diverges fast enough that the score is essentially linear on ℐmI_m, recovering the classical m m behavior. Theorem 5.3. Assume Assumption 5.2 and fix α∈(0,1]α∈(0,1]. With the tuning above and with X~m X_m defined as the unique central root of g^m(x)=0 g_m(x)=0 in ℐmI_m, there exists an event ℰmE_m with ℙ(ℰm)→1P(E_m)→ 1 such that, on ℰmE_m, g^m g_m has a unique zero in ℐmI_m; moreover, X~m→μ X_m pμ, and m(X~m−μ)⇒N(0,σ2) m( X_m-μ) N(0,σ^2). Figure 4: Comparison of the effectiveness of the three-stage interaction model,The grey-filled area represents the experimental results without the three-stage interaction, while the orange part shows the results obtained with the three-stage interaction. It can be seen that the three-stage interactive solution is effective. Figure 5: The performance of the proposed method under different noise conditions is illustrated: the left figure shows contamination by noise groups of varying proportions, while the right figure depicts contamination ratios within the same group. It can be observed that the proposed method demonstrates strong noise resistance. Figure 6: Training reward curves of ARE and FlightGPT. 5.2 Heavy-Tailed Regime Assumption 5.4 ((1+ϵ)(1+ε)-moment). Fix ϵ∈(0,1]ε∈(0,1]. Assume [X1]=μE[X_1]=μ and [|X1−μ|1+ϵ]≤v1+ϵ<∞E[ |X_1-μ |^1+ε]≤ v_1+ε<∞. For any confidence level δ∈(0,1/2)δ∈(0,1/2), choose cm(δ):=τ(v1+ϵmlog⁡(2/δ))11+ϵ,τ>0,c_m(δ):=τ ( v_1+ε\,m (2/δ) ) 11+ε,τ>0, and take any rm→∞r_m→∞ with rm=o(cm(δ))r_m=o(c_m(δ)) (e.g., rm=cm(δ)1/2r_m=c_m(δ)^1/2). Here cm(δ)c_m(δ) is tuned to balance (controlled) bias from downweighting large residuals with concentration of the empirical score. Theorem 5.5. Assume Assumption 5.4 and fix α∈(0,1]α∈(0,1]. With cm=cm(δ)c_m=c_m(δ) and the central-root estimator X~m X_m, there exist constants C1,C2>0C_1,C_2>0 depending only on α,ϵ,τα,ε,τ such that for all sufficiently large m, ℙ(|X~m−μ|≤C1(v1+ϵlog⁡(C2/δ)m)ϵ1+ϵ)≥1−δ.P ( | X_m-μ |≤ C_1 ( v_1+ε (C_2/δ)m ) ε1+ε )≥ 1-δ. This single-block deviation bound is the basic ingredient for median-of-means aggregation, which amplifies constant-success probability to an arbitrary confidence level. 5.3 Asymptotic Theory for ARE We now aggregate across blocks using a sample median. In the finite-variance regime, a block-level CLT combined with the classical median CLT yields a n n limit and the usual π/2π/2 variance inflation, giving a constant asymptotic relative efficiency. Theorem 5.6. Assume Assumption 5.2 and fix α∈(0,1]α∈(0,1]. Let k=kn→∞k=k_n→∞ and kn=o(n)k_n=o(n), and set m:=⌊n/kn⌋→∞m:= n/k_n →∞. On each block, compute X~j X_j using the finite-variance tuning in Theorem 5.3, namely for a fixed γ>0γ>0, cm:=m1/2+γ,rm:=mγ/4.c_m:=m^1/2+γ,r_m:=m^γ/4. Assume in addition that the block-level distribution of X~j X_j admits the local normal approximation needed for the sample-median CLT (formalized in Appendix G). Then n(μ~−μ)⇒N(0,π2σ2). n\,( μ-μ)\ \ N\! (0, π2σ^2 ). Consequently, relative to the sample mean X¯n X_n, the asymptotic relative efficiency is 2π. 2π. Theorem 5.7 (Heavy-tailed deviation bound). Assume Assumption 5.4 with some ϵ∈(0,1]ε∈(0,1] and fix α∈(0,1]α∈(0,1]. Let δ∈(0,1/2)δ∈(0,1/2) and choose the number of blocks k:=⌈8log⁡(2δ)⌉,m:=⌊n/k⌋.k:= 8 \! ( 2δ ) ,m:= n/k . On each block, run the central-root estimator X~j X_j with the scale cm:=cm(δ0)=τ(v1+ϵmlog⁡(2/δ0))11+ϵ,δ0:=14,c_m:=c_m( _0)=τ ( v_1+ε\,m (2/ _0) ) 11+ε, _0:= 14, and any radius rm→∞r_m→∞ such that rm=o(cm)r_m=o(c_m) (e.g., rm=cm1/2r_m=c_m^1/2). Then there exist constants C1,C2>0C_1,C_2>0 depending only on α,ϵ,τα,ε,τ such that for all sufficiently large n, ℙ(|μ~−μ|≤C1(v1+ϵlog⁡(C2/δ)n)ϵ1+ϵ)≥1−δ.P\! ( | μ-μ |≤ C_1 ( v_1+ε (C_2/δ)n ) ε1+ε )≥ 1-δ. 6 Experiments 6.1 Mathematical Reasoning Experiments In this section, we evaluate the proposed ARE and DACR on challenging mathematical reasoning benchmarks across LLMs of different parameter scales and model families. We train on MATH 7.5k (Hendrycks et al., 2021) and test on MATH 500 (Lightman et al., 2023) (in-distribution), using AMC23, Gaokao2023en (Zhang et al., 2024), and MinervaMath (Lewkowycz et al., 2022) as out-of-distribution suites to assess generalization. We compare against a baseline that matches our collaboration setting but replaces ARE/DACR with standard GRPO. We organize the empirical study around three research questions, examining overall effectiveness under noisy rewards, the contribution of the three-stage interaction, and robustness under diverse contamination patterns. Q1: Is the proposed multi-agent robust training effective? To evaluate robustness under noisy rewards, we inject heavy-tailed, peaked Cauchy noise into the reward model and train on mathematical reasoning datasets. As shown in Table 1, our method consistently improves accuracy over the GRPO-based baseline on both ID and OOD tests. The gains can be substantial; e.g., on AMC23, Qwen2.5-1.5B-Instruct improves from 22.5%22.5\% to 32.5%32.5\% (+10%10\%), suggesting the proposed multi-agent robust training is effective across model scales and families. Q2: Is the proposed three-stage interaction effective? We ablate the interaction protocol by training two agents that answer independently while keeping all other settings unchanged. Results in Figure 4 show that the three-stage interaction yields consistently better reasoning performance across LLMs; in particular, on AMC23, Qwen2.5-7B-Instruct improves from 47.5%47.5\% to 60%60\% (+12.5%12.5\%). Q3: How robust is the method to noise? We further stress-test robustness under diverse contamination patterns (Figure 5): (i) group-level corruption where a fraction of groups are selected and 20% of samples within each selected group are contaminated, and (i) within-group corruption where the contamination ratio varies inside the same group. Both settings evaluate accuracy on MATH500. Across all conditions, our method degrades gracefully and remains consistently more robust to reward contamination than the baseline. Table 1: Comparison of accuracy rates of various methods on the mathematical reasoning dataset Model Method MATH500 AIME24 AIME25 AMC23 Gaokao 2023en Minerva Math Qwen2.5 (1.5B-Instruct) Base Model 35.80 - - 30.00 23.50 5.150 baseline-GroupA 56.00 - - 22.50 45.50 13.20 baseline-GroupC 55.00 3.300 - 22.50 45.50 13.60 ARE-GroupA 56.60 3.300 - 32.50 46.00 14.00 ARE-GroupC 55.60 3.300 3.300 27.50 45.50 15.40 Llama3.1 (8B-Instruct) Base Model 32.60 - - 5.000 3.390 7.350 baseline-GroupB 49.40 3.300 - 20.00 39.70 17.60 baseline-GroupC 49.20 6.700 3.300 22.50 39.20 17.30 ARE-GroupB 51.00 6.700 - 22.50 42.60 19.90 ARE-GroupC 50.40 6.700 - 27.50 40.30 18.80 Qwen2.5 (7B-Instruct) Base Model 68.00 6.667 10.00 37.50 42.04 18.38 baseline-GroupA 77.60 16.70 6.700 50.00 62.60 27.60 baseline-GroupB 77.40 13.30 10.00 52.50 61.30 28.70 ARE-GroupA 77.60 10.00 10.00 60.00 61.80 28.30 ARE-GroupB 77.40 13.30 10.00 57.70 61.60 30.50 6.2 Embodied VLN Evaluation Figure 7: Visualization of the navigation task. The UAV trajectory, predicted target, and ground-truth target are annotated. To assess whether ARE transfers beyond pure language reasoning, we evaluated it on an embodied aerial VLN task by applying ARE to a visual language model. We built on FlightGPT (Cai et al., 2025) and replaced its standard GRPO optimization with ARE, fine-tuning a Qwen2-VL-2B-Instruct backbone (Wang et al., 2024) with downscaled image resolutions to fit available compute; due to current constraints, we focused on validating the Adaptive Robust Estimator component, leaving full multi-agent interaction to future work. In VLN, the agent follows a natural-language instruction D and visual observations from a 3D environment E, starting from (0,θ0)(p_0, _0) and executing discrete actions at∈a_t ; an episode succeeds if the final position Tp_T is within a threshold δ of the target. Following (Cai et al., 2025), we reported four standard metrics: Navigation Error (NE), the Euclidean distance between the final position and the ground-truth target (lower is better); Success Rate (SR), the fraction of episodes whose final position is within 20 meters of the target (higher is better); Oracle Success Rate (OSR), the fraction of episodes that reach within 20 meters of the target at any time (higher is better); and Success weighted by Path Length (SPL), which penalizes unnecessarily long trajectories (higher is better). Figure 6 shows that ARE achieved consistently higher training rewards than FlightGPT, suggesting improved optimization stability under the same VLN setup; although we used a smaller backbone than in (Cai et al., 2025) (leading to lower absolute performance), the comparison to the GRPO baseline remained controlled. Consistently, Table 2 shows that ARE outperformed both the base model and FlightGPT across all splits, with the largest gains on unseen environments (e.g., NE 145.56!→! 124.38 on valunseenval_unseen and 140.78!→! 129.71 on testunseentest_unseen) while also improving SR/OSR/SPL, indicating better localization, higher success, and more efficient trajectories. Table 2: Aerial VLN results (higher is better for SR/OSR/SPL; lower is better for NE). Split Method NE ↓ SR ↑ OSR ↑ SPL ↑ VS Base 193.92 0.3300 1.0000 0.0000 FlightGPT 148.19 2.6700 13.670 0.0200 ARE 139.49 3.0000 14.330 0.0300 VUS Base 190.84 1.3300 1.3300 0.0100 FlightGPT 145.56 3.6700 18.000 0.0300 ARE 124.38 6.3300 18.330 0.0600 TUS Base 195.95 1.3300 2.3300 0.0100 FlightGPT 140.78 2.3300 17.330 0.0200 ARE 129.71 4.6700 18.330 0.0400 • VS: val_seen, VUS: val_unseen, TUS: test_unseen. 7 Conclusion This work addressed a key instability in MARL for LLM reasoning: noisy, heavy-tailed rewards can make GRPO’s batch-mean advantage normalization fragile, leading to oscillation or divergence. We proposed a robust collaborative training framework that coupled a structured DACR protocol with a cross-improvement reward for more reliable credit assignment, and introduced ARE to replace the batch mean in GRPO-style advantage estimation. We established theoretical guarantees for ARE under both finite-variance and heavy-tailed regimes, and experiments on mathematical reasoning benchmarks showed consistent gains in accuracy and training stability across model scales and families, including homogeneous and heterogeneous agent pairs, especially under heavy-tailed reward noise. We further validated ARE in an embodied aerial VLN setting by replacing FlightGPT’s GRPO with ARE on a Qwen2-VL-2B VLM, yielding higher training rewards and improved NE/SR/OSR/SPL across all splits, with the largest gains on unseen environments. Impact Statement This paper presents work whose goal is to advance the field of machine learning by developing robust multi-agent policy optimization methods for collaborative reasoning under noisy, heavy-tailed rewards. We do not feel any potential societal consequences of this work necessary to be discussed here. References Barron (2019) Barron, J. T. A general and adaptive robust loss function. CVPR, 2019. Black & Anandan (1996) Black, M. J. and Anandan, P. The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields. Computer vision and image understanding, 63(1):75–104, 1996. Black & Rangarajan (1996) Black, M. J. and Rangarajan, A. On the unification of line processes, outlier rejection, and robust statistics with applications in early vision. Int. J. Comput. Vision, 19(1):57–91, July 1996. ISSN 0920-5691. Cai et al. (2025) Cai, H., Dong, J., Tan, J., Deng, J., Li, S., Gao, Z., Wang, H., Su, Z., Sumalee, A., and Zhong, R. Flightgpt: Towards generalizable and interpretable uav vision-and-language navigation with vision-language models, 2025. URL https://arxiv.org/abs/2505.12835. Chebrolu et al. (2021) Chebrolu, N., Läbe, T., Vysotska, O., Behley, J., and Stachniss, C. Adaptive robust kernels for non-linear least squares problems. IEEE Robotics and Automation Letters, 6(2):2240–2247, 2021. Chen et al. (2025a) Chen, Z., Ai, T., Li, Y., Li, G., Wei, Y., Zhou, W., Li, G., Yu, B., Chen, Z., Sun, H., Zhuang, F., Li, J., Wang, D., and Ban, Y. Llmboost: Make large language models stronger with boosting, 2025a. URL https://arxiv.org/abs/2512.22309. Chen et al. (2025b) Chen, Z., Ji, Z., Mao, Q., Wu, H., Cheng, J., Qin, B., Li, Z., Li, J., Sun, K., Wang, Z., et al. Scoring, reasoning, and selecting the best! ensembling large language models via a peer-review process. arXiv preprint arXiv:2512.23213, 2025b. Chen et al. (2025c) Chen, Z., Li, J., Chen, P., Li, Z., Sun, K., Luo, Y., Mao, Q., Li, M., Xiao, L., Yang, D., et al. Harnessing multiple large language models: A survey on llm ensemble. arXiv preprint arXiv:2502.18036, 2025c. Chen et al. (2026) Chen, Z., Li, G., Ai, T., Li, Y., Huang, Z., Zhou, W., Zhuang, F., Liu, X., Li, J., Wang, D., et al. Weak-driven learning: How weak agents make strong agents stronger. arXiv preprint arXiv:2602.08222, 2026. Dennis Jr & Welsch (1978) Dennis Jr, J. E. and Welsch, R. E. Techniques for nonlinear least squares and robust regression. Communications in Statistics-simulation and Computation, 7(4):345–359, 1978. Geman & Geman (1986) Geman, D. and Geman, S. Bayesian image analysis. In Disordered systems and biological organization, p. 301–319. Springer, 1986. Hastie et al. (2001) Hastie, T., Tibshirani, R., and Friedman, J. The Elements of Statistical Learning. Springer Series in Statistics. Springer New York Inc., New York, NY, USA, 2001. He et al. (2025) He, X., Ban, Y., Zou, J., Wei, T., Cook, C., and He, J. Llm-forest: Ensemble learning of llms with graph-augmented prompts for data imputation. In Findings of the Association for Computational Linguistics: ACL 2025, p. 6921–6936, 2025. Hendrycks et al. (2021) Hendrycks, D., Burns, C., Kadavath, S., Arora, A., Basart, S., Tang, E., Song, D., and Steinhardt, J. Measuring mathematical problem solving with the math dataset, 2021. URL https://arxiv.org/abs/2103.03874. Hitchcox & Forbes (2022) Hitchcox, T. and Forbes, J. R. Mind the gap: Norm-aware adaptive robust loss for multivariate least-squares problems. IEEE Robotics and Automation Letters, 7(3):7116–7123, 2022. Hu (2025) Hu, J. Reinforce++: A simple and efficient approach for aligning large language models. arXiv preprint arXiv:2501.03262, 2025. Huang et al. (2026a) Huang, Z., Xia, X., Ren, Y., Zheng, J., Wang, X., Zhang, Z., Xie, H., Liang, S., Chen, Z., Xiao, X., et al. Does your reasoning model implicitly know when to stop thinking? arXiv preprint arXiv:2602.08354, 2026a. Huang et al. (2026b) Huang, Z., Xia, X., Ren, Y., Zheng, J., Xiao, X., Xie, H., Huaqiu, L., Liang, S., Dai, Z., Zhuang, F., et al. Real-time aligned reward model beyond semantics. arXiv preprint arXiv:2601.22664, 2026b. Huber (2011) Huber, P. J. Robust statistics. In International encyclopedia of statistical science, p. 1248–1251. Springer, 2011. Humbert et al. (2022) Humbert, P., Le Bars, B., and Minvielle, L. Robust kernel density estimation with median-of-means principle. In International Conference on Machine Learning, p. 9444–9465. PMLR, 2022. Jung et al. (2024) Jung, K., Hitchcox, T., and Forbes, J. R. An adaptive graduated nonconvexity loss function for robust nonlinear least-squares solutions. IEEE Transactions on Robotics, 2024. Leclerc (1989) Leclerc, Y. G. Constructing simple stable descriptions for image partitioning. International journal of computer vision, 3(1):73–102, 1989. Lecué et al. (2020) Lecué, G., Lerasle, M., and Mathieu, T. Robust classification via mom minimization. Machine learning, 109(8):1635–1665, 2020. Lewkowycz et al. (2022) Lewkowycz, A., Andreassen, A., Dohan, D., Dyer, E., Michalewski, H., Ramasesh, V., Slone, A., Anil, C., Schlag, I., Gutman-Solo, T., Wu, Y., Neyshabur, B., Gur-Ari, G., and Misra, V. Solving quantitative reasoning problems with language models, 2022. URL https://arxiv.org/abs/2206.14858. Liao et al. (2025) Liao, J., Wen, M., Wang, J., and Zhang, W. Marft: Multi-agent reinforcement fine-tuning. arXiv preprint arXiv:2504.16129, 2025. Lightman et al. (2023) Lightman, H., Kosaraju, V., Burda, Y., Edwards, H., Baker, B., Lee, T., Leike, J., Schulman, J., Sutskever, I., and Cobbe, K. Let’s verify step by step, 2023. URL https://arxiv.org/abs/2305.20050. Lin et al. (2025) Lin, H., Cao, S., Wang, S., Wu, H., Li, M., Yang, L., Zheng, J., and Qin, C. Interactive learning for llm reasoning. arXiv preprint arXiv:2509.26306, 2025. Liu et al. (2025a) Liu, H., Fu, Z., Ding, M., Ning, R., Zhang, C., Liu, X., and Zhang, Y. Logical reasoning in large language models: A survey. arXiv preprint arXiv:2502.09100, 2025a. Liu et al. (2025b) Liu, S., Chen, T., Liang, Z., Lyu, X., and Amato, C. Llm collaboration with multi-agent reinforcement learning. arXiv preprint arXiv:2508.04652, 2025b. Liu et al. (2025c) Liu, W., Hu, H., Zhou, J., Ding, Y., Li, J., Zeng, J., He, M., Chen, Q., Jiang, B., Zhou, A., et al. Mathematical language models: A survey. ACM Computing Surveys, 58(6):1–37, 2025c. Lu et al. (2022) Lu, L., Yu, Y., de Lamare, R. C., and Yang, X. Tukey’s biweight m-estimate with conjugate gradient adaptive learning. IEEE Signal Processing Letters, 29:1117–1121, 2022. Lu et al. (2026) Lu, X., Wang, X., Chai, J., Yin, G., Lin, W., Chen, Z., Luo, Y., Zhuang, F., Ban, Y., and Wang, D. Contextual rollout bandits for reinforcement learning with verifiable rewards. arXiv preprint arXiv:2602.08499, 2026. Ma et al. (2024) Ma, H., Hu, T., Pu, Z., Boyin, L., Ai, X., Liang, Y., and Chen, M. Coevolving with the other you: Fine-tuning llm with sequential cooperative multi-agent reinforcement learning. Advances in Neural Information Processing Systems, 37:15497–15525, 2024. Motwani et al. (2024) Motwani, S. R., Smith, C., Das, R. J., Rafailov, R., Laptev, I., Torr, P. H., Pizzati, F., Clark, R., and de Witt, C. S. Malt: Improving reasoning with multi-agent llm training. arXiv preprint arXiv:2412.01928, 2024. Park et al. (2025) Park, C., Han, S., Guo, X., Ozdaglar, A. E., Zhang, K., and Kim, J.-K. Maporl: Multi-agent post-co-training for collaborative large language models with reinforcement learning. In Proceedings of the 63rd Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), p. 30215–30248, 2025. Rafailov et al. (2023) Rafailov, R., Sharma, A., Mitchell, E., Manning, C. D., Ermon, S., and Finn, C. Direct preference optimization: Your language model is secretly a reward model. Advances in neural information processing systems, 36:53728–53741, 2023. Ren et al. (2025) Ren, T., Jiang, J., Yang, H., Tian, W., and Peng, Y. RiskPO: Risk-based policy optimization with verifiable reward for LLM post-training. In NeurIPS 2025 Workshop MLxOR: Mathematical Foundations and Operational Integration of Machine Learning for Uncertainty-Aware Decision-Making, 2025. URL https://openreview.net/forum?id=8hxqmh25ZH. Rieder (2012) Rieder, H. Robust Asymptotic Statistics: Volume I. Springer Science & Business Media, 2012. Shao et al. (2024) Shao, Z., Wang, P., Zhu, Q., Xu, R., Song, J., Bi, X., Zhang, H., Zhang, M., Li, Y. K., Wu, Y., and Guo, D. Deepseekmath: Pushing the limits of mathematical reasoning in open language models, 2024. URL https://arxiv.org/abs/2402.03300. Shao et al. (2025) Shao, Z., Luo, Y., Lu, C., Ren, Z., Hu, J., Ye, T., Gou, Z., Ma, S., and Zhang, X. Deepseekmath-v2: Towards self-verifiable mathematical reasoning. arXiv preprint arXiv:2511.22570, 2025. Sun (2021) Sun, Q. Do we need to estimate the variance in robust mean estimation? arXiv preprint arXiv:2107.00118, 2021. Wan et al. (2025) Wan, Z., Li, Y., Wen, X., Song, Y., Wang, H., Yang, L., Schmidt, M., Wang, J., Zhang, W., Hu, S., et al. Rema: Learning to meta-think for llms with multi-agent reinforcement learning. arXiv preprint arXiv:2503.09501, 2025. Wang et al. (2024) Wang, P., Bai, S., Tan, S., Wang, S., Fan, Z., Bai, J., Chen, K., Liu, X., Wang, J., Ge, W., Fan, Y., Dang, K., Du, M., Ren, X., Men, R., Liu, D., Zhou, C., Zhou, J., and Lin, J. Qwen2-vl: Enhancing vision-language model’s perception of the world at any resolution. arXiv preprint arXiv:2409.12191, 2024. Wang et al. (2025) Wang, P.-Y., Liu, T.-S., Wang, C., Li, Z., Wang, Y., Yan, S., Jia, C., Liu, X.-H., Chen, X., Xu, J., et al. A survey on large language models for mathematical reasoning. ACM Computing Surveys, 2025. Wang & Xu (2025) Wang, X. and Xu, M. Multi-agent multi-armed bandit with fully heavy-tailed dynamics. arXiv preprint arXiv:2501.19239, 2025. Xue et al. (2025) Xue, X., Zhou, Y., Zhang, G., Zhang, Z., Li, Y., Zhang, C., Yin, Z., Torr, P., Ouyang, W., and Bai, L. Comas: Co-evolving multi-agent systems via interaction rewards. arXiv preprint arXiv:2510.08529, 2025. Yang et al. (2026) Yang, F., Chen, Z., Wang, X., Lu, X., Chai, J., Yin, G., Lin, W., Ma, S., Zhuang, F., Wang, D., et al. Your group-relative advantage is biased. arXiv preprint arXiv:2601.08521, 2026. Zhang et al. (2024) Zhang, X., Li, C., Zong, Y., Ying, Z., He, L., and Qiu, X. Evaluating the performance of large language models on gaokao benchmark, 2024. URL https://arxiv.org/abs/2305.12474. Zhang et al. (2026) Zhang, Z., Huang, Z., Xia, X., Wang, D., Zhuang, F., Ma, S., Ding, N., Yang, Y., Li, J., and Ban, Y. Heterogeneous agent collaborative reinforcement learning. arXiv preprint arXiv:2603.02604, 2026. Zou et al. (2025) Zou, J., Ban, Y., Li, Z., Qi, Y., Qiu, R., Yang, L., and He, J. Transformer copilot: Learning from the mistake log in LLM fine-tuning. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. URL https://openreview.net/forum?id=MRvxlTlkNQ. Appendix A Related Work A.1 Multi-Agent Collaboration and Training Recent advances in large language models have inspired substantial interest in multi-agent collaboration(Chen et al., 2025a; He et al., 2025; Chen et al., 2026, 2025b; Lu et al., 2026), where multiple LLM-based agents interact to solve tasks that exceed the capabilities of a single model. Early explorations primarily focused on role specialization, task decomposition, and collaborative protocols such as planning, debating, and mutual verification. These frameworks improve solution quality through interaction but typically rely on fixed, frozen models without joint optimization. More recent studies extend this paradigm by treating multi-agent systems as interactive learning environments, in which agents coordinate, compete, or exchange information while being jointly optimized through reinforcement learning. Collectively, these efforts demonstrate that coordinated multi-agent LLMs can exhibit stronger reasoning, more robust generalization, and higher reliability than single-agent systems. A representative example is CORY (Ma et al., 2024), which formulates two symbiotic agents: a pioneer and an observer, both initialized from the same base model. The pioneer generates an initial response, while the observer synthesizes an improved refinement based on the pioneer’s output. The two agents periodically switch roles to foster cooperative evolution. MARFT (Liao et al., 2025) models multi-agent collaboration as a sequential decision-making process that aligns with the PPO training paradigm. This approach allows agents to share a common backbone while adopting distinct capabilities through different LoRA adapters. Beyond these architectures, several works explicitly incorporate reinforcement learning for multi-agent reasoning. Lin et al. (2025) propose an interactive learning framework where agents dynamically choose cooperative or competitive behaviors based on problem difficulty and model capability. Their three-stage discussion protocol enhances information exchange, and the reward design incorporates distributional properties of other agents’ returns. Training is performed using GRPO, ultimately improving the reasoning ability of individual agents through multi-agent synergy. Similarly, CoMAS (Xue et al., 2025) constructs a rich interaction loop based on solution proposal, peer evaluation, and scoring, and employs an LLM-as-a-judge mechanism to convert discussion signals into reward feedback. These rewards are then optimized via the REINFORCE++ (Hu, 2025) algorithm, enabling decentralized and scalable emergent learning without external supervision.(Yang et al., 2026) provides a principled theoretical analysis of group-based advantage estimation. Several recent works extend GRPO into explicit multi-agent regimes. ReMA (Wan et al., 2025) introduces meta-level reasoning into hierarchical multi-agent training, involving high-level and low-level agents trained under both shared and independent GRPO setups with single- and multi-round interaction schemes. MAGRPO (Liu et al., 2025b) further generalizes GRPO to collaborative fine-tuning across multiple agents, where a centralized group-relative advantage guides joint optimization while retaining decentralized execution; this method has demonstrated effectiveness in program synthesis tasks. MAPoRL (Park et al., 2025) proposes a post-training framework that combines multi-agent reinforcement learning with a carefully designed reward mechanism to explicitly train LLMs for coordinated decision-making in complex environments. Meanwhile, MALT (Motwani et al., 2024) constructs a division-of-labor multi-agent system with iterative collaboration and experience-based improvement, leveraging tree-structured data generation, value-iteration-style credit assignment, and a hybrid SFT–DPO (Rafailov et al., 2023) post-training pipeline to overcome the limitations of single-agent reasoning. While these works have established a solid foundation for multi-agent reinforcement learning with LLMs, they rarely address the robustness challenges of GRPO in multi-agent optimization. In contrast, our approach introduces a new collaboration mechanism and explicitly incorporates robust reinforcement learning techniques to stabilize multi-agent GRPO training. A.2 Robust Estimation Methods Robust estimation methods are designed to maintain statistical reliability when data deviate from idealized assumptions, such as in the presence of outliers, model misspecification, or heavy-tailed noise (Huber, 2011; Rieder, 2012). Classical estimators—e.g., the sample mean or ordinary least squares—can be highly sensitive to a small fraction of corrupted observations, leading to unstable inference and degraded predictive performance. After decades of development, robust estimation methods can be categorized—according to their evolutionary trajectory—into the design of robust loss functions, adaptive procedures, tuning-free approaches, and median-of-means (MOM) estimators. The core idea of robust loss functions is to minimize the influence of abnormally large per-sample residuals or gradients on the overall objective (Jung et al., 2024). Representative choices include the Huber (Huber, 2011), pseudo-Huber (Hastie et al., 2001), Cauchy (Black & Anandan, 1996), Tukey’s biweight (Lu et al., 2022), and Geman–McClure (GM) (Geman & Geman, 1986) losses . A key caveat, however, is the risk of over-robustness: overly aggressive down-weighting of extreme observations may suppress rare yet genuine signals. Moreover, the landscape of robust losses is broad, and selecting an appropriate form typically relies heavily on expert knowledge and empirical tuning. To mitigate this issue, Barron (2019) proposed an adaptive loss (as shown in (9)) which, via a shape parameter α governing tail behavior and robustness and a scale parameter c controlling the width of the near-zero quadratic basin, unifies many classical losses within a single continuously adjustable—and potentially automatically learnable—framework. Nevertheless, both classes of approaches may suffer from substantial computational overhead for tuning-parameter selection (e.g., common choices such as Lepski’s method and cross-validation can be computationally expensive), and the resulting objectives are often nonconvex (Jung et al., 2024), which can lead to poor local optima, sensitivity to initialization, and training instability. To address these two major challenges, a third line of work—tuning-free methods—has emerged. The key idea is to determine the value of the tuning parameter (typically depending on the target estimation accuracy and the population second moment of the data) under which the resulting estimator can achieve a sub-Gaussian rate of convergence. A representative example is the work of Sun (2021), which builds on the pseudo-Huber loss and derives an explicit criterion for choosing the tuning parameter. A fourth class of approaches that has attracted considerable attention in recent years is the MoM estimator (Lecué et al., 2020; Humbert et al., 2022). The procedure first partitions the sample into multiple groups, computes the sample mean within each group, and then takes the median of these group-wise means. Under the mere existence of a finite second moment, the MoM estimator can achieve sub-Gaussian rates of convergence. In Appendix B, we further discuss that, under both random contamination and adversarial contamination, the MoM estimator may still fail to be robust; moreover, attaining strong empirical performance typically requires a careful trade-off between the number of blocks and the sample size within each block. Appendix B Discussion on the MOM Estimator and its Robustness In this section, we review the Median-of-Means (MoM) estimator as our baseline and discuss its key limitations under two contamination regimes: random and adversarial contamination. Let X1,…,XnX_1,…,X_n be i.i.d. real-valued random variables with mean μ=[Xi]μ=E[X_i]. The empirical mean X¯=1n∑i=1nXi X= 1n _i=1^nX_i is statistically optimal under light-tailed distributions (e.g., sub-Gaussian), yet it can be highly unstable in the presence of heavy tails or contaminated samples. As noted in Section 1, MoM provides a canonical robust alternative by combining blockwise averaging with a robust aggregation step. Concretely, fix k≥2k≥ 2 and randomly partition [n][n] into k disjoint blocks B1,…,BkB_1,…,B_k of size m=⌊n/k⌋m= n/k . Define the block means X¯j=1|Bj|∑i∈BjXi,j=1,…,k, X_j= 1|B_j| _i∈ B_jX_i, j=1,…,k, and aggregate them via the median: μ^=median⁡(X¯1,…,X¯k). μ=median( X_1,…, X_k). This construction exhibits a transparent robustness mechanism: blockwise averaging reduces the variance contributed by the uncontaminated observations, while the median discards the influence of a minority of corrupted blocks. In particular, under the mild assumption σ2=Var(Xi)<∞σ^2=Var(X_i)<∞, MoM achieves a sub-Gaussian-type concentration bound. For any δ∈(0,1)δ∈(0,1), choosing k≃c0log⁡(1/δ)k c_0 (1/δ) yields ℙ(|μ^−μ|≤C0σlog⁡(1/δ)n)≥1−δ,P\! ( | μ-μ |≤ C_0σ (1/δ)n )≥ 1-δ, for universal constants c0,C0>0c_0,C_0>0 (Sun, 2021). Hence, even without sub-Gaussian tails, MoM retains the optimal n−1/2n^-1/2 rate and the canonical log⁡(1/δ) (1/δ) confidence scaling. Practically, k controls the robustness–efficiency trade-off: increasing k improves tolerance to corruption through the median at the cost of smaller block sizes. A common default is k=⌈c0log⁡(1/δ)⌉k= c_0 (1/δ) (with c0c_0 in the range 22–88); when δ is not specified, the heuristic k≈⌈log⁡n⌉k≈ n is often used. We now turn to a simple analysis and discussion. In random contamination, assume each observation is independently corrupted with probability κ∈(0,1)κ∈(0,1). Randomly partition the n samples into k disjoint blocks, each of size m=⌊n/k⌋.m= n/k . Call a block clean if none of its m samples is contaminated. Under independent contamination, the probability that a fixed block is clean is p0=Pr⁡(a block is clean)=(1−κ)m,p_0= (a block is clean)=(1-κ)^m, where we ignore the minor imbalance due to the floor operation and possible leftovers. Let Xj=block j is cleanX_j=1\block $j$ is clean\ for j=1,…,kj=1,…,k. Then Xj∼Bernoulli(p0)X_j (p_0) i.i.d., and the number of clean blocks is S:=∑j=1kXj∼Binomial(k,p0).S:= _j=1^kX_j (k,p_0). A key mechanism behind MoM is that the median is taken over a majority of blocks; thus, if more than half of the blocks are clean (i.e., S>k/2S>k/2), then the median is computed from predominantly uncontaminated block means. This motivates the condition p0>1/2p_0>1/2. When p0>1/2p_0>1/2, a Chernoff–Hoeffding (KL-form) bound yields Pr⁡(S≤k/2)≤exp⁡(−kKL(1/2∥p0)), (S≤ k/2)≤ \! (-k\,KL (1/2\,\|\,p_0 ) ), where the binary KL divergence is KL(1/2∥p0)=12ln⁡1/2p0+12ln⁡1/21−p0.KL (1/2\,\|\,p_0 )= 12 \! 1/2p_0+ 12 \! 1/21-p_0. Hence, provided p0>1/2p_0>1/2, the failure event “at most half of the blocks are clean” has probability exponentially small in k. The condition p0>1/2p_0>1/2 is equivalent to (1−κ)m>12⟺m<ln⁡2−ln⁡(1−κ).(1-κ)^m> 12 m< 2- (1-κ). For small κ, using −ln⁡(1−κ)≈κ- (1-κ)≈κ, this becomes the convenient rule of thumb m≲ln⁡2κ.m 2κ. Thus, smaller blocks (smaller m, larger k) make it easier to ensure that a majority of blocks are clean, which is precisely the regime in which MoM can retain sub-Gaussian-type deviation bounds under only finite second moments. It is important to note that “smaller blocks are cleaner” does not mean “the smaller the better.” Decreasing m increases p0=(1−κ)mp_0=(1-κ)^m and improves the probability of having a clean majority, but it also increases the variance of each within-block mean because fewer samples are averaged within a block. Therefore, obtaining good estimation accuracy requires balancing the number of blocks k and the block size m. In the adversarial contamination, an attacker can strategically place corrupted observations across blocks. In particular, it suffices to inject a single extremely large outlier into each of ⌈k/2⌉ k/2 distinct blocks to arbitrarily shift the MoM output. Indeed, within any contaminated block of size m, the block mean takes the form X¯block=1m∑i=1mXi, X_block= 1m _i=1^mX_i, so a single outlier of magnitude M shifts the block mean by Δ=Mm. = Mm. Although the multiplicative factor 1/m1/m attenuates the effect, the attacker may choose M arbitrarily large, hence Δ can be made arbitrarily large. Consequently, once at least half of the blocks are contaminated, the median over block means can be driven to an arbitrary value. This yields a finite-sample breakdown point on the order of ϵbd≈⌈k/2⌉n, _bd≈ k/2 n, since contaminating ⌈k/2⌉ k/2 blocks requires only ⌈k/2⌉ k/2 adversarial points among n observations. In particular, if k is fixed (or grows sublinearly in n), then ϵbd→0 _bd→ 0 as n→∞n→∞, indicating extremely poor robustness in the strict breakdown-point sense. Equivalently, even when the overall contamination proportion is very small, MoM can be completely compromised if more than half of the blocks each contain at least one outlier. Finally, the choice of k exhibits an intrinsic tension. If k is too large, then m=⌊n/k⌋m= n/k becomes too small and the variance of each block mean inflates; if k is too small, then blocks are too large and it becomes easier for contamination to hit a majority of blocks (especially under random contamination). Therefore, MoM entails a structural trade-off between k and m that is difficult to calibrate in practice. Appendix C Simulation Studies of the ARE In the main text, we validate the effectiveness of the ARE estimator through several concrete LLM applications. In this section, we further demonstrate its performance via simulation studies. We consider one-dimensional robust mean estimation with ground-truth mean μ⋆=1μ =1. For each sample size n∈200,500,1000,2000,5000n∈\200,500,1000,2000,5000\, we generate i.i.d. samples X1,…,XnX_1,…,X_n under five data-generating processes (DGP1–DGP5) that progressively increase distributional difficulty from light tails to heavy tails and explicit corruption. DGP1 (Gaussian) serves as a light-tailed baseline: Xi=μ⋆+σZi,Zi∼(0,1),X_i=μ +σ Z_i, Z_i (0,1), with default σ=1σ=1. DGP2 (centered log-normal) captures skewed heavy tails while preserving the mean: letting Yi∼(0,1)Y_i (0,1), Xi=μ⋆+(exp⁡(τYi)−exp⁡(τ2/2)),X_i=μ + ( (τ Y_i)- (τ^2/2) ), where τ>0τ>0 controls tail-heaviness and skewness (default τ=1.5τ=1.5). DGP3 (Student-t) models symmetric heavy tails with finite variance: Xi=μ⋆+sTi,Ti∼tν,ν>2,X_i=μ +sT_i, T_i t_ν, ν>2, with default (s,ν)=(1,4)(s,ν)=(1,4); in our implementation Ti=Zi/Vi/νT_i=Z_i/ V_i/ν with Zi∼(0,1)Z_i (0,1) and Vi∼χ2(ν)V_i χ^2(ν). DGP4 (symmetric Pareto) probes the infinite-variance regime while keeping the mean finite: Xi=μ⋆+SiUi,X_i=μ +S_iU_i, where Si∈−1,+1S_i∈\-1,+1\ is Rademacher (equiprobable) and Ui≥1U_i≥ 1 has Pareto tail ℙ(Ui>u)=u−αP(U_i>u)=u^-α for u≥1u≥ 1, with α∈(1,2)α∈(1,2) (default α=1.5α=1.5). Finally, DGP5 (contamination) explicitly injects outliers on top of a clean baseline Xi(0)∼(μ⋆,σ02)X_i^(0) (μ , _0^2) with default σ0=1 _0=1. Under random (mean-preserving) contamination, independently for each i we set Xi=Xi(0)X_i=X_i^(0) with probability 1−κ1-κ, and with probability κ replace Xi←μ⋆+SiM,X_i←μ +S_iM, where Si∈−1,+1S_i∈\-1,+1\ is equiprobable (default (κ,M)=(0.05,100)(κ,M)=(0.05,100)). Under adversarial one-sided contamination, after drawing the clean sample we replace exactly ⌊κn⌋ κ n points by the outlier μ⋆+Mμ +M. To further stress-test median-of-means (MoM) style aggregation, we also consider a block-aware adversary: we deterministically partition the first mkmk samples into k blocks of equal size m=⌊n/k⌋m= n/k (discarding leftovers), then select r=⌈k/2⌉r= k/2 blocks and inject one outlier μ⋆+Mμ +M into each selected block (optionally capped by a global budget ⌊κn⌋ κ n ). We compare the estimators l2, Cauchy, GM, Adapt, AMB (adaptive Maxwell–Boltzmann (Hitchcox & Forbes, 2022)), their GNC variants GNC_Adapt, GNC_AMB, GNC_TLS (truncated least squares), and the MoM-integrated variant ARE. Each method is initialized at the sample mean μ(0)=1n∑i=1nXiμ^(0)= 1n _i=1^nX_i and returns a scalar estimate μ μ. For each pair (DGP,n)(DGP,n), we run R=100R=100 Monte Carlo replications. To ensure fair comparisons, we use common random numbers: in replication r∈0,…,R−1r∈\0,…,R-1\ we set the seed to base_seed+r base\_seed+r (default base_seed=12345 base\_seed=12345), generate one dataset, and evaluate all estimators on the same data. Performance is summarized by the error e=μ^−μ⋆e= μ-μ through MAE and MSE, MAE=[|e|],MSE=[e2],MAE=E[|e|], =E[e^2], together with robust quantiles of |e||e| (Median, Q90, and Q95 computed via linear interpolation on sorted absolute errors). Tables 3–7 report the experimental results for DGP1–DGP5 across all sample sizes and estimators. Table 3: Results on DGP1 (Gaussian): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 0.057287 0.074298 0.082726 0.089879 0.086309 0.089576 0.086304 0.067855 0.093291 MSE 0.005173 0.008747 0.010770 0.012688 0.011817 0.012600 0.011710 0.007302 0.013051 Median 0.047519 0.062271 0.068645 0.074465 0.070920 0.074508 0.072363 0.056955 0.080019 Q90 0.117608 0.152255 0.173105 0.185062 0.181983 0.185578 0.178264 0.140570 0.187690 Q95 0.137931 0.199400 0.216909 0.229500 0.224137 0.225030 0.220215 0.163779 0.224390 500 MAE 0.034494 0.044633 0.049594 0.053918 0.052523 0.053909 0.052433 0.041655 0.063132 MSE 0.001838 0.002937 0.003645 0.004326 0.004069 0.004322 0.004044 0.002647 0.006196 Median 0.029021 0.040605 0.043591 0.046434 0.045928 0.046208 0.046618 0.035880 0.052503 Q90 0.069328 0.089981 0.100150 0.106283 0.105652 0.106068 0.105136 0.085231 0.132819 Q95 0.081448 0.099868 0.114402 0.126156 0.119561 0.126863 0.119259 0.097621 0.152423 1000 MAE 0.024141 0.032439 0.036162 0.039273 0.038195 0.039322 0.038274 0.029353 0.044525 MSE 0.000877 0.001590 0.001982 0.002343 0.002228 0.002344 0.002231 0.001330 0.003116 Median 0.020658 0.028489 0.031086 0.034594 0.034696 0.034421 0.034927 0.025598 0.037707 Q90 0.048351 0.064046 0.073582 0.079218 0.079125 0.079847 0.079148 0.059341 0.089249 Q95 0.057840 0.077597 0.085396 0.092410 0.090472 0.092422 0.089938 0.068942 0.108791 2000 MAE 0.018659 0.025467 0.028122 0.030252 0.029671 0.030232 0.029631 0.022884 0.035059 MSE 0.000534 0.000973 0.001192 0.001395 0.001326 0.001394 0.001323 0.000812 0.001860 Median 0.016846 0.022130 0.024437 0.026106 0.026217 0.026044 0.026388 0.020445 0.029204 Q90 0.038356 0.052066 0.056480 0.061172 0.058984 0.061170 0.058737 0.046448 0.072581 Q95 0.044255 0.059942 0.066329 0.071238 0.069291 0.071240 0.068512 0.055781 0.082386 5000 MAE 0.011756 0.015152 0.016491 0.017618 0.017328 0.017611 0.017323 0.014490 0.020988 MSE 0.000215 0.000365 0.000439 0.000507 0.000490 0.000507 0.000489 0.000317 0.000683 Median 0.010020 0.012336 0.013943 0.014923 0.014729 0.014922 0.014766 0.013010 0.017922 Q90 0.023354 0.032167 0.035496 0.037974 0.037047 0.037991 0.036445 0.027374 0.041582 Q95 0.027684 0.036030 0.040371 0.043431 0.043515 0.043434 0.042303 0.033452 0.049732 Table 4: Results on DGP2 (Centered LogNormal): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 0.440732 2.383597 2.428866 2.437117 2.441007 2.457052 2.348575 2.228842 2.442037 MSE 0.319861 5.684898 5.905568 5.977074 5.964692 6.039979 5.522237 4.972983 5.968493 Median 0.374759 2.385970 2.433266 2.458548 2.446208 2.460386 2.348119 2.225001 2.446614 Q90 0.835407 2.453506 2.501665 2.523363 2.513925 2.521825 2.452893 2.326485 2.530520 Q95 1.092987 2.474307 2.514922 2.537888 2.530098 2.537393 2.473458 2.343185 2.551792 500 MAE 0.284348 2.384501 2.431850 2.456983 2.446077 2.457022 2.342375 2.230460 2.456446 MSE 0.139719 5.686893 5.914873 6.037721 5.984414 6.037914 5.489974 4.976731 6.035749 Median 0.233371 2.385287 2.431715 2.457309 2.446462 2.457377 2.335733 2.230171 2.458230 Q90 0.599346 2.426119 2.470482 2.496481 2.488803 2.496624 2.425090 2.284673 2.504184 Q95 0.708262 2.434107 2.482159 2.508471 2.499906 2.508461 2.444517 2.296450 2.517809 1000 MAE 0.211270 2.385299 2.432640 2.457747 2.448701 2.457779 2.329159 2.231928 2.457179 MSE 0.077294 5.690247 5.918278 6.041038 5.996701 6.041196 5.426657 4.982442 6.038572 Median 0.172539 2.384034 2.432140 2.457033 2.448279 2.457028 2.324848 2.231803 2.457498 Q90 0.424567 2.417506 2.463839 2.489640 2.481998 2.489626 2.387608 2.270945 2.494162 Q95 0.530385 2.423898 2.470188 2.493565 2.486350 2.493547 2.406430 2.283688 2.506638 2000 MAE 0.156645 2.385716 2.433062 2.458169 2.449900 2.458201 2.318107 2.231655 2.458427 MSE 0.043052 5.691973 5.920085 6.042870 6.002305 6.043028 5.374561 4.980825 6.044297 Median 0.121270 2.386073 2.433150 2.457977 2.449910 2.458020 2.313820 2.231244 2.459061 Q90 0.320521 2.409282 2.456259 2.480855 2.473205 2.480893 2.358750 2.261899 2.485424 Q95 0.394085 2.417107 2.461496 2.485910 2.478026 2.485956 2.375311 2.270079 2.491608 5000 MAE 0.105697 2.385195 2.432413 2.457464 2.449503 2.457496 2.307456 2.231414 2.458300 MSE 0.018055 5.689290 5.916752 6.039241 6.000178 6.039394 5.324652 4.979428 6.043399 Median 0.091034 2.385240 2.432084 2.457331 2.449294 2.457379 2.306014 2.231975 2.458453 Q90 0.216724 2.399676 2.445936 2.470685 2.462998 2.470644 2.330204 2.251057 2.474467 Q95 0.262409 2.405306 2.450735 2.475080 2.466920 2.475065 2.338020 2.255622 2.477549 Table 5: Results on DGP3 (Student-t): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 0.079329 0.077027 0.085139 0.091804 0.089159 0.091806 0.088340 0.075760 0.106217 MSE 0.010140 0.009076 0.011105 0.013011 0.012191 0.013010 0.011964 0.009004 0.017967 Median 0.065749 0.066426 0.074241 0.079843 0.077080 0.079817 0.076304 0.065225 0.086817 Q90 0.165889 0.156022 0.175618 0.185670 0.181116 0.185723 0.181239 0.155254 0.216114 Q95 0.197331 0.181987 0.201629 0.219082 0.212402 0.219154 0.211227 0.181042 0.265095 500 MAE 0.053192 0.050686 0.055188 0.059119 0.057425 0.059175 0.057249 0.050691 0.070144 MSE 0.004389 0.003959 0.004725 0.005423 0.005150 0.005426 0.005107 0.004040 0.007499 Median 0.047555 0.042452 0.048275 0.049438 0.048677 0.049481 0.048193 0.041320 0.063075 Q90 0.108891 0.100059 0.111604 0.121137 0.118760 0.120968 0.118494 0.103940 0.144723 Q95 0.124714 0.122820 0.130709 0.143039 0.141928 0.143011 0.136721 0.123777 0.162166 1000 MAE 0.036055 0.034389 0.036687 0.038872 0.038216 0.038801 0.038277 0.035842 0.045915 MSE 0.002049 0.001818 0.002119 0.002405 0.002328 0.002397 0.002328 0.001978 0.003375 Median 0.030219 0.029908 0.031209 0.032017 0.031530 0.031926 0.031413 0.030449 0.039262 Q90 0.076757 0.071994 0.077547 0.081623 0.079673 0.081599 0.080059 0.072346 0.096297 Q95 0.087859 0.085480 0.095375 0.100148 0.098976 0.100240 0.098355 0.083406 0.117685 2000 MAE 0.027005 0.025412 0.027685 0.029742 0.029301 0.029659 0.029244 0.026076 0.035955 MSE 0.001130 0.001019 0.001203 0.001371 0.001331 0.001361 0.001327 0.001082 0.002001 Median 0.024092 0.021646 0.023765 0.025259 0.025824 0.025264 0.025472 0.022653 0.029967 Q90 0.053873 0.050953 0.057300 0.062042 0.060372 0.061152 0.059646 0.054660 0.073035 Q95 0.063653 0.059921 0.065569 0.071133 0.069951 0.070878 0.070029 0.063390 0.083720 5000 MAE 0.016767 0.015550 0.016721 0.017673 0.017497 0.017613 0.017462 0.016209 0.022485 MSE 0.000438 0.000380 0.000446 0.000505 0.000494 0.000503 0.000493 0.000417 0.000775 Median 0.014579 0.013322 0.014540 0.015093 0.015065 0.015025 0.015005 0.013681 0.019660 Q90 0.033488 0.031232 0.033738 0.035248 0.035180 0.035245 0.035151 0.033712 0.045893 Q95 0.041146 0.037998 0.041529 0.044344 0.043462 0.044339 0.043444 0.041594 0.051865 Table 6: Results on DGP4 (Symmetric Pareto): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 0.412362 1.438328 1.518497 1.522176 1.517095 1.500790 1.564941 1.595852 1.326587 MSE 0.590463 2.072604 2.413169 2.459166 2.413946 2.255401 2.458244 2.645237 1.767638 Median 0.253832 1.436787 1.502406 1.498419 1.502898 1.501094 1.564170 1.650789 1.325492 Q90 0.848063 1.514769 1.572911 1.565906 1.574786 1.566497 1.664394 1.763932 1.438845 Q95 1.155875 1.537037 1.602260 1.596185 1.600452 1.599216 1.694309 1.795882 1.461990 500 MAE 0.346482 1.439051 1.505764 1.513167 1.516621 1.497524 1.574814 1.638352 1.339161 MSE 0.441919 2.072407 2.268632 2.356217 2.371738 2.243820 2.482337 2.718181 1.798361 Median 0.214237 1.435626 1.503844 1.501487 1.503984 1.495642 1.576356 1.655792 1.345156 Q90 0.679670 1.488627 1.554075 1.548618 1.551994 1.545145 1.633025 1.723029 1.425720 Q95 1.001716 1.505555 1.565078 1.559460 1.565327 1.556368 1.648180 1.741259 1.447595 1000 MAE 0.404736 1.435188 1.568055 1.499761 1.502030 1.496583 1.581408 1.650191 1.172465 MSE 6.035140 2.060631 4.486813 2.249927 2.260321 2.244539 2.502444 2.742182 1.592772 Median 0.174446 1.433928 1.502667 1.497949 1.503413 1.497246 1.585403 1.658140 1.350682 Q90 0.617388 1.473958 1.538876 1.532762 1.537795 1.533496 1.626974 1.709116 1.419330 Q95 0.977243 1.485230 1.548406 1.543379 1.549137 1.545707 1.641047 1.721744 1.435046 2000 MAE 0.277339 1.483650 1.553332 1.549032 1.553667 1.498893 1.583741 1.640998 1.358529 MSE 1.570783 3.414067 3.633623 3.623838 3.619150 2.247003 2.509354 2.721683 1.850885 Median 0.134661 1.434549 1.503553 1.499047 1.504036 1.497680 1.589020 1.660958 1.364235 Q90 0.499724 1.461581 1.527771 1.521979 1.528136 1.521673 1.621502 1.693832 1.421070 Q95 0.723037 1.468161 1.534138 1.527576 1.533661 1.527547 1.628375 1.702626 1.432571 5000 MAE 0.172371 1.433381 1.503061 1.498617 1.504033 1.498603 1.585251 1.610063 1.376493 MSE 0.291315 2.054726 2.259312 2.245967 2.262231 2.245924 2.518299 2.668343 1.896178 Median 0.095233 1.433980 1.502981 1.498358 1.503860 1.498899 1.592539 1.660890 1.377615 Q90 0.312437 1.447404 1.516828 1.511865 1.517150 1.511360 1.608947 1.678992 1.423776 Q95 0.446499 1.452532 1.519658 1.515645 1.521086 1.515568 1.613800 1.682779 1.436759 Table 7: Results on DGP5 (Random contamination): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 1.275762 0.076195 0.085042 0.092851 0.088960 0.092902 0.082093 0.069203 0.101623 MSE 2.589767 0.009235 0.011385 0.013490 0.012538 0.013503 0.010758 0.007540 0.016209 Median 1.034393 0.063265 0.071296 0.077451 0.074361 0.077476 0.064493 0.058816 0.085344 Q90 2.562352 0.154705 0.171824 0.187716 0.182085 0.187708 0.170956 0.143705 0.219979 Q95 3.040650 0.198719 0.218241 0.232770 0.224837 0.232767 0.207915 0.179176 0.251584 500 MAE 0.803141 0.045911 0.051108 0.055784 0.054351 0.055812 0.051871 0.043391 0.065892 MSE 0.999150 0.003058 0.003778 0.004492 0.004241 0.004495 0.003863 0.002878 0.006678 Median 0.647632 0.042344 0.047195 0.050444 0.050155 0.050443 0.046949 0.037627 0.055327 Q90 1.648696 0.088117 0.097799 0.107246 0.102712 0.107400 0.100482 0.087705 0.133173 Q95 1.987630 0.103240 0.113667 0.125033 0.119585 0.124933 0.113398 0.103239 0.162270 1000 MAE 0.559301 0.033795 0.037628 0.040745 0.039726 0.040750 0.038179 0.030827 0.047724 MSE 0.478545 0.001704 0.002122 0.002517 0.002390 0.002517 0.002215 0.001423 0.003510 Median 0.493729 0.028837 0.032411 0.035934 0.034792 0.035898 0.033473 0.027246 0.038480 Q90 1.097783 0.068452 0.077436 0.084156 0.082611 0.084147 0.076178 0.059305 0.097886 Q95 1.313985 0.080211 0.088760 0.097425 0.095241 0.097411 0.092335 0.073964 0.115841 2000 MAE 0.402444 0.025673 0.028418 0.030721 0.030015 0.030728 0.029719 0.023478 0.036463 MSE 0.252472 0.000999 0.001222 0.001433 0.001362 0.001434 0.001334 0.000846 0.002092 Median 0.342546 0.022232 0.024743 0.026295 0.026008 0.026305 0.026590 0.020136 0.030277 Q90 0.814902 0.052501 0.057734 0.063119 0.060448 0.063057 0.059638 0.046563 0.077859 Q95 0.983756 0.063525 0.067954 0.074476 0.071649 0.074469 0.071677 0.054432 0.091296 5000 MAE 0.254088 0.015730 0.017206 0.018488 0.018124 0.018487 0.018104 0.014815 0.023329 MSE 0.099986 0.000389 0.000470 0.000546 0.000526 0.000546 0.000520 0.000330 0.000835 Median 0.214689 0.013191 0.014152 0.015257 0.014829 0.015312 0.014979 0.013506 0.020551 Q90 0.528963 0.031760 0.035809 0.038966 0.038035 0.038986 0.038560 0.027583 0.047545 Q95 0.641245 0.037724 0.041986 0.045497 0.044011 0.045536 0.043942 0.032985 0.055688 Table 8: Results on DGP5 (Adversarial contamination): robust error summaries over R=100R=100 replications. n Metric l2 Cauchy GM Adapt AMB GNC_Adapt GNC_AMB GNC_TLS ARE 200 MAE 4.998228 0.076277 0.085139 0.093145 0.089440 0.093207 0.072658 0.069398 0.103908 MSE 24.987188 0.009264 0.011434 0.013583 0.012628 0.013600 0.008495 0.007742 0.016552 Median 4.998399 0.063285 0.072447 0.080556 0.077372 0.080661 0.059623 0.057831 0.091392 Q90 5.084834 0.158220 0.175776 0.190045 0.185091 0.190049 0.154349 0.143999 0.212980 Q95 5.108672 0.194569 0.210418 0.232594 0.216446 0.232717 0.181021 0.173504 0.255916 500 MAE 4.999549 0.045390 0.050404 0.054917 0.053448 0.054955 0.044563 0.042151 0.067477 MSE 24.997235 0.003090 0.003836 0.004572 0.004302 0.004578 0.003019 0.002751 0.007088 Median 4.999719 0.040502 0.044703 0.046358 0.046212 0.046366 0.037025 0.036002 0.055582 Q90 5.050312 0.089395 0.100666 0.109910 0.105790 0.110204 0.091927 0.086210 0.142707 Q95 5.064844 0.105286 0.118548 0.130789 0.127427 0.130905 0.106573 0.102324 0.159420 1000 MAE 5.001430 0.033238 0.037159 0.040485 0.039282 0.040500 0.032143 0.030005 0.047662 MSE 25.015113 0.001671 0.002089 0.002487 0.002350 0.002489 0.001573 0.001385 0.003576 Median 5.001633 0.028801 0.031912 0.034680 0.034012 0.034688 0.028010 0.027005 0.039794 Q90 5.037464 0.067114 0.075948 0.083690 0.080102 0.083690 0.064081 0.060579 0.096793 Q95 5.044215 0.077899 0.088521 0.096025 0.093886 0.096265 0.075479 0.071472 0.119724 2000 MAE 5.001119 0.025654 0.028436 0.030799 0.030083 0.030811 0.024935 0.023116 0.036394 MSE 25.011676 0.001013 0.001251 0.001476 0.001393 0.001478 0.000961 0.000838 0.002057 Median 5.002635 0.020993 0.023428 0.025935 0.025570 0.025954 0.021426 0.019531 0.032832 Q90 5.029097 0.052707 0.058848 0.063274 0.060947 0.063288 0.052107 0.050658 0.073951 Q95 5.036158 0.063923 0.069545 0.075759 0.074667 0.075907 0.059192 0.057579 0.086389 5000 MAE 5.001120 0.015621 0.016972 0.018272 0.017900 0.018279 0.015694 0.014902 0.021548 MSE 25.011411 0.000387 0.000465 0.000539 0.000520 0.000539 0.000384 0.000336 0.000751 Median 5.001860 0.013595 0.014565 0.015719 0.015379 0.015783 0.013764 0.013436 0.017448 Q90 5.019530 0.031161 0.034034 0.036961 0.035849 0.036985 0.030413 0.028134 0.044297 Q95 5.022983 0.037216 0.040165 0.044437 0.043214 0.044466 0.037985 0.035315 0.054321 Across the five data-generating processes (DGP1–DGP5), the dominant pattern is that ARE offers its clearest benefits precisely in the regimes where robustness is genuinely needed: skewed heavy tails, infinite-variance tails, and explicit contamination. In the skewed heavy-tailed DGP2 (centered log-normal) and the infinite-variance DGP4 (symmetric Pareto with tail index ξ∈(1,2)ξ∈(1,2)), several non-MoM robust losses (e.g., Cauchy/GM/Adapt/AMB and their GNC counterparts) exhibit a pronounced non-vanishing error floor: their MAE and Median remain essentially flat as n grows, with high quantiles (Q90/Q95) also failing to improve materially. In contrast, ARE consistently suppresses the tail risk in these settings, yielding markedly smaller Q90/Q95 and, in DGP4, substantially improved overall error summaries at representative sample sizes (e.g., a striking reduction in MAE/MSE at n=1000n=1000). This behavior is consistent with the two-layer robustness of the method: the inner GNC mechanism stabilizes the optimization of a nonconvex robust objective, while the outer median-of-means aggregation truncates the influence of a minority of blocks containing extreme observations, thereby preventing the estimator from being dominated by rare but catastrophic samples. The advantage becomes even more transparent under explicit contamination (DGP5). Under random contamination, the sample mean (l2) suffers a dramatic degradation in all metrics, while ARE maintains small errors and noticeably better tail quantiles, indicating strong protection against sporadic outliers. Under adversarial contamination, l2 is effectively “locked” by the injected outliers: its MAE stays near 55 and does not meaningfully decrease with n, reflecting a complete breakdown of classical averaging. In sharp contrast, ARE remains accurate and continues to improve with n, with particularly strong control of Q90/Q95, which highlights its robustness to worst-case corruption consistent with MoM-type guarantees (so long as a majority of blocks remain uncontaminated). Importantly, in the benign baselines DGP1 (Gaussian) and DGP3 (Student-t with finite variance), ARE incurs only a modest constant-factor efficiency loss relative to the best light-tail methods (e.g., l2 or GNC_TLS) and does not exhibit instability. Overall, these results suggest that ARE trades a small amount of efficiency in easy regimes for substantial and reliable gains in heavy-tailed and contaminated regimes, with its most consistent improvements appearing in the high-quantile error metrics (Q90/Q95) that directly capture tail-risk robustness. Appendix D Introduction to Adaptive Losses. In this paper, we employ the adaptive loss function introduced by Barron (2019), defined as follows: ρ(ϵ;α,c)=|α−2|α((1+(ϵ/c)2|α−2|)α/2−1),ρ(ε;α,c)= |α-2|α ( (1+ (ε/c)^2|α-2| )^α/2-1 ), (9) where α∈(−∞,2]α∈(-∞,2] is a shape parameter controlling the loss’s robustness, while c>0c>0 is a scale parameter governing the quadratic bowl near ϵ=0ε=0. Through simple analysis, we can immediately recover several common losses—such as the L2, Charbonnier/pseudo-Huber/L1-L2, Cauchy/Lorentzian (Black & Anandan, 1996), Geman–McClure (Geman & Geman, 1986), and Welsch/Leclerc losses (Dennis Jr & Welsch, 1978; Leclerc, 1989)—as special cases of the adaptive loss ρ(ϵ;α,c)ρ(ε;α,c) in (9) by setting α to particular values, as summarized as follows, ρ(ϵ;α,c)=12(ϵ/c)2α=2(ϵ/c)2+1−1α=1log⁡(12(ϵ/c)2+1)α=02(ϵ/c)2/((ϵ/c)2+4)α=−21−exp⁡(−12(ϵ/c)2)α=−∞,∂ρ∂ϵ(ϵ;α,c)=ϵ/c2α=2(ϵ/c2)/(ϵ/c)2+1α=12ϵ/(ϵ2+2c2)α=0(ϵ/c2)((ϵ/c)2/4+1)−3/2α=−2(ϵ/c2)exp⁡(−12(ϵ/c)2)α=−∞ρ(ε;α,c)= cases 12(ε/c)^2&α=2\\ (ε/c)^2+1-1&α=1\\ ( 12(ε/c)^2+1 )&α=0\\ 2(ε/c)^2/ ((ε/c)^2+4 )&α=-2\\ 1- (- 12(ε/c)^2 )&α=-∞ cases, ∂ρ∂ε(ε;α,c)= casesε/c^2&α=2\\ (ε/c^2)/ (ε/c)^2+1&α=1\\ 2ε/(ε^2+2c^2)&α=0\\ (ε/c^2) ((ε/c)^2/4+1 )^-3/2&α=-2\\ (ε/c^2) (- 12(ε/c)^2 )&α=-∞ cases Figure 8: Adaptive robust loss ρ(ϵ;α,c)ρ(ε;α,c) and its gradient under different parameter settings. From left to right and top to bottom, the four panels show: (top-left) the loss ρ(ϵ;α,c)ρ(ε;α,c) for varying shape α with fixed c=1c=1; (top-right) the corresponding gradient ∂ρ/∂ϵ∂ρ/∂ε for varying α with fixed c=1c=1; (bottom-left) the loss for varying scale c with fixed α=0α=0; (bottom-right) the corresponding gradient for varying c with fixed α=0α=0. The adaptive loss ρ(ϵ;α,c)ρ(ε;α,c) provides a unified two-parameter family that spans a broad spectrum of robust penalties through a shape parameter α∈(−∞,2]α∈(-∞,2] and a scale parameter c>0c>0. As illustrated in Figure 8, decreasing α continuously interpolates from the quadratic ℓ2 _2 loss (α=2α=2), which assigns unbounded influence to large residuals, to increasingly robust losses whose growth slows down for large |ϵ||ε|. This continuum makes α a direct knob for robustness: smaller α reduces the marginal penalty on large errors and thus mitigates the effect of outliers or heavy-tailed noise. The scale parameter c controls the width of the near-zero quadratic regime. For fixed α (e.g., α=0α=0), larger c enlarges the approximately quadratic region around ϵ=0ε=0, thereby behaving more like least squares on moderate residuals and improving statistical efficiency under mild noise. Conversely, smaller c triggers earlier deviation from quadratic growth, increasing robustness by treating moderate-to-large residuals more conservatively. Hence, c mediates a local efficiency–robustness trade-off by determining the transition point at which the loss begins to downweight large errors. The gradients ∂ρ/∂ϵ∂ρ/∂ε make the robustness mechanism explicit through the induced influence function. When α is small, the gradient decays rapidly as |ϵ||ε| grows and can approach zero for very large residuals, effectively suppressing the contribution of extreme observations during optimization. This bounded (or rapidly vanishing) influence stabilizes gradient-based updates and prevents optimization from being dominated by a small number of corrupted samples. The scale c further adjusts the decay rate of ∂ρ/∂ϵ∂ρ/∂ε: smaller c yields faster saturation (stronger outlier rejection but potentially slower progress when many residuals are moderately large), whereas larger c preserves larger gradients over a wider range (often faster convergence locally but weaker robustness). Overall, the (α,c)(α,c) parameterization enables data-dependent control of robustness versus efficiency within a single coherent framework, which is particularly useful for robust estimation and learning under heavy-tailed or contaminated data. Appendix E The GNC Paradigm for Mean Estimation in Subsection 4.2 This appendix details the graduated nonconvexity (GNC) procedure used to solve the nonconvex subproblem (5) (equivalently (8)) after obtaining (α(1),c(1))(α^(1),c^(1)). We present (i) a GNC surrogate for the adaptive loss, (i) the induced IRLS weights and the associated Black–Rangarajan (BR) regularization, and (i) a complete GNC-IRLS algorithm with practical continuation schedules and stopping criteria. Let f(β,α(1))f(β,α^(1)) be a shape mapping controlled by a continuation parameter β. We construct a one-parameter family of surrogate losses by replacing the shape parameter α(1)α^(1) with f(β,α(1))f(β,α^(1)): ρβ(ϵ;α(1),c(1))≔ρ(ϵ;f(β,α(1)),c(1)), _β(ε;α^(1),c^(1))\; \;ρ\! (ε;\ f(β,α^(1)),\ c^(1) ), (10) where ρ(⋅;α,c)ρ(·;α,c) denotes the adaptive robust loss used in the main text. To avoid ambiguity and ensure numerical stability, we recommend writing ρβ _β in the standard Barron form with a scale parameter: ρβ(ϵ;α(1),c(1))=|f−2|f[((ϵ/c(1))2|f−2|+1)f2−1],f≡f(β,α(1)). _β(ε;α^(1),c^(1))= |f-2|f [ ( (ε/c^(1))^2|f-2|+1 ) f2-1 ], f≡ f(β,α^(1)). (11) This definition guarantees the correct quadratic limit as f→2f→ 2 and matches the commonly used scaling convention (cf. Barron-style adaptive losses). If a different scaling convention is adopted in the main text, (11) should be adjusted accordingly; the derivations below remain identical up to constant factors. A valid shape mapping must satisfy limβ→β1f(β,α(1))=2,limβ→β2f(β,α(1))=α(1),f(β,α(1))≤2, _β→ _1f(β,α^(1))=2, _β→ _2f(β,α^(1))=α^(1), f(β,α^(1))≤ 2, (12) where β1,β2 _1, _2 are surrogate-dependent constants (or endpoints) determined by the chosen parameterization of β. Intuitively, f(β,α(1))→2f(β,α^(1))→ 2 yields a convex quadratic surrogate, while f(β,α(1))→α(1)f(β,α^(1))→α^(1) recovers the original adaptive loss (typically nonconvex). We consider three concrete choices of f(β,α(1))f(β,α^(1)): Example 1 (Polynomial decay; β:∞→1β:\ ∞→ 1). f1(β,α(1))=2−2−α(1)βp,p>0,f_1(β,α^(1))=2- 2-α^(1)β^p, p>0, (13) so that β→∞β→∞ implies f1→2f_1→ 2 and β→1β→ 1 implies f1→α(1)f_1→α^(1). Example 2 (Exponential transition; β: 0+→∞β:\ 0^+→∞). f2(β,α(1))=α(1)e−1/β+2(1−e−1/β),f_2(β,α^(1))=α^(1)e^-1/β+2 (1-e^-1/β ), (14) so that β→0+β→ 0^+ implies f2→2f_2→ 2 and β→∞β→∞ implies f2→α(1)f_2→α^(1). Example 3 (Rational form; β: 0+→∞β:\ 0^+→∞). f3(β,α(1))=2+α(1)βq1+βq,q>0,f_3(β,α^(1))= 2+α^(1)β^q1+β^q, q>0, (15) so that β→0+β→ 0^+ implies f3→2f_3→ 2 and β→∞β→∞ implies f3→α(1)f_3→α^(1). Figure 9: Three instantiations of the interpolation function f(β,α(1))f(β,α^(1)): (a) polynomial decay f1=2−2−α(1)βpf_1=2- 2-α^(1)β^p with β:∞→1β:∞→ 1; (b) exponential transition f2=α(1)e−1/β+2(1−e−1/β)f_2=α^(1)e^-1/β+2(1-e^-1/β) with β:0+→∞β:0^+→∞; (c) rational form f3=2+α(1)βq1+βqf_3= 2+α^(1)β^q1+β^q with β:0+→∞β:0^+→∞, all smoothly interpolating between the quadratic regime (f=2f=2) and the adaptive regime (f=α(1)f=α^(1)) for different values of α(1)∈−2,0,1α^(1)∈\-2,0,1\. Let ϵi=Xi−X _i=X_i-X and consider the surrogate objective minX∈ℝ⁡1m∑i=1mρβ(ϵi;α(1),c(1)). _X \ 1m _i=1^m _β( _i;α^(1),c^(1)). (16) At fixed (β,α(1),c(1))(β,α^(1),c^(1)), IRLS solves a sequence of weighted least-squares problems minX∈ℝ⁡12m∑i=1mwiϵi2, _X \ 12m _i=1^mw_i\, _i^2, (17) with weights obtained by matching derivatives w.r.t. X: wi=1ϵi∂ϵiρβ(ϵi;α(1),c(1)).w_i\;=\; 1 _i ∂ _i _β( _i;α^(1),c^(1)). (18) With ρβ _β defined in (11), straightforward differentiation yields wi=1(c(1))2((ϵi/c(1))2|f−2|+1)f2−1,f=f(β,α(1)).w_i= 1(c^(1))^2 ( ( _i/c^(1))^2|f-2|+1 ) f2-1, f=f(β,α^(1)). (19) For numerical robustness, it is often convenient to use the following special-case expressions when f is close to singular values (and otherwise use (19)): wi=1(c(1))2,f=2,2ϵi2+2(c(1))2,f=0,1(c(1))2((ϵi/c(1))2|f−2|+1)f2−1,otherwise,w_i= cases 1(c^(1))^2,&f=2,\\[4.0pt] 2 _i^2+2(c^(1))^2,&f=0,\\[6.0pt] 1(c^(1))^2 ( ( _i/c^(1))^2|f-2|+1 ) f2-1,&otherwise, cases (20) and to clip wiw_i to (0,1](0,1] when required by the BR-duality constraint. Optimizing (17) jointly over (X,(wi)i=1m)(X,(w_i)^m_i=1) without regularization degenerates to wi=0w_i=0 for all i. To prevent this, we exploit the Black–Rangarajan duality (Black & Rangarajan, 1996). Under standard conditions (e.g., concavity of ϕ(z)≔ρβ(c(1)z)φ(z)\! \! _β(c^(1) z) in z), there exists an outlier-process penalty Φρβ(⋅) _ _β(·) such that ρβ(ϵ)=minw∈(0,1]⁡12wϵ2+Φρβ(w). _β(ε)= _w∈(0,1] 12wε^2+ _ _β(w). (21) Plugging (21) into (16) yields the BR-regularized formulation minX∈ℝ,wi∈(0,1]⁡12m∑i=1m(wiϵi2+Φρβ(wi)), _X ,\,w_i∈(0,1]\ 12m _i=1^m (w_i _i^2+ _ _β(w_i) ), (22) which matches (8) with the surrogate loss ρβ _β. At fixed X, minimizing (22) over wiw_i recovers the IRLS update (18). For fixed weights, (17) admits a closed-form minimizer: X←∑i=1mwiXi∑i=1mwi.X← _i=1^mw_iX_i _i=1^mw_i. (23) The GNC continuation progressively changes β so that f(β,α(1))f(β,α^(1)) moves from the quadratic regime (f≈2f≈ 2) to the target regime (f≈α(1)f≈α^(1)). Example 1 schedule (β:∞→1β:∞→ 1). We initialize β0 _0 at a large value (e.g., β0≫1 _0 1 so that f1(β0,α(1))≈2f_1( _0,α^(1))≈ 2) and decrease it geometrically toward 11: βk=1+βk−1−1γ,γ>1. _k=1+ _k-1-1γ, γ>1. (24) This update ensures βk↓1 _k 1 and hence f1(βk,α(1))↓α(1)f_1( _k,α^(1)) α^(1). Example 2/3 schedule (β:0+→∞β:0^+→∞). We initialize β0 _0 small (e.g., β0≪1 _0 1 so that f2/3(β0,α(1))≈2f_2/3( _0,α^(1))≈ 2) and increase it geometrically: βk=γβk−1,γ>1, _k=γ _k-1, γ>1, (25) so that βk↑∞ _k ∞ and f2/3(βk,α(1))→α(1)f_2/3( _k,α^(1))→α^(1). We terminate the outer continuation once |f(β,α(1))−α(1)|≤εf, |f(β,α^(1))-α^(1) |≤ _f, (26) or after a maximum number of continuation steps. Algorithm 1 summarizes the full procedure. The method alternates IRLS updates of (X,(wi)i=1m)(X,(w_i)^m_i=1) at a fixed β, and then updates β to increase nonconvexity. This yields a robust solver with improved basins of attraction compared to directly optimizing the original nonconvex loss. We use geometric schedules (24) (Example 1) or (25) (Examples 2–3) with γ∈[1.2,2]γ∈[1.2,2] in our experiments. For initialization, we set X to the sample mean (or median for additional robustness) and choose β0 _0 such that f(β0,α(1))f( _0,α^(1)) is sufficiently close to 22 (e.g., |f(β0,α(1))−2|≤10−2|f( _0,α^(1))-2|≤ 10^-2). Each inner iteration costs (m)O(m) due to the closed-form updates (19) and (23). Hence the total cost is (mTmaxKmax)O(m\,T_ \,K_ ). The continuation parameter β controls the rate at which nonconvexity is injected. A faster schedule (larger γ in (24)–(25)) can accelerate convergence but may reduce robustness to poor initialization; conversely, a slower schedule yields more stable iterations with larger attraction basins at the expense of additional outer steps. Algorithm 1 GNC-IRLS for solving (8) (Step-1 update) 0: Measurements (Xi)i=1m(X_i)_i=1^m; parameters (α(1),c(1))(α^(1),c^(1)); mapping f(β,α(1))f(β,α^(1)); schedule β←(β)β (β); tolerances εX,εf _X, _f; iteration limits TmaxT_ (inner) and KmaxK_ (outer). 0: X~(1) X^(1). 1: X←1m∑i=1mXiX← 1m _i=1^mX_i (or median); choose β←β0β← _0 such that f(β0,α(1))≈2f( _0,α^(1))≈ 2. 2: for k=1,2,…,Kmaxk=1,2,…,K_ do 3: f←f(β,α(1))f← f(β,α^(1)). 4: for t=1,2,…,Tmaxt=1,2,…,T_ do 5: ϵi←Xi−X _i← X_i-X, wi←1ϵi∂ϵiρ(ϵi;f,c(1)),w_i← 1 _i ∂ _iρ( _i;f,c^(1)), wi←min⁡1,max⁡δ,wi,∀iw_i← \1, \δ,w_i\\,\ ∀ i 6: Xnew←∑i=1mwiXi∑i=1mwiX_new← _i=1^mw_iX_i _i=1^mw_i 7: if |Xnew−X|≤εX|X_new-X|≤ _X then 8: X←XnewX← X_new; 9: end if 10: X←XnewX← X_new. 11: end for 12: if |f(β,α(1))−α(1)|≤εf |f(β,α^(1))-α^(1) |≤ _f then 13: break 14: end if 15: β←(β)β (β) (increase nonconvexity). 16: end for 17: return X~(1)←X X^(1)← X. Appendix F Technical Proofs F.1 Proof of Lemma 4.1 Proof. We divide the analysis into two cases. When α>2α>2, we have |α−2|=α−2|α-2|=α-2, so ρ(ϵ;α,c)=α−2α[(1+(ϵ/c)2α−2)α/2−1].ρ(ε;α,c)= α-2α [ (1+ (ε/c)^2α-2 )^α/2-1 ]. Let y=1+(ϵ/c)2α−2(≥1).y=1+ (ε/c)^2α-2 (≥ 1). Differentiating with respect to α (via the product and chain rules) yields ∂ρ∂α=2α2(yα/2−1)+α−22αyα/2log⁡y−12yα/2−1(y−1). ∂ρ∂α= 2α^2 (y^α/2-1 )+ α-22α\,y^α/2 y- 12\,y^α/2-1(y-1). For fixed α, ϵε, and c, define the function F:[1,∞)→ℝF:[1,∞) by F(y):=∂ρ∂α(ϵ,α,c)=2α2(yα/2−1)+α−22αyα/2log⁡y−12yα/2−1(y−1).F(y):= ∂ρ∂α(ε,α,c)= 2α^2 (y^α/2-1 )+ α-22α\,y^α/2 y- 12\,y^α/2-1(y-1). Direct differentiation in y gives dyF(y)=α−24yα/2−2(1−y+ylog⁡y)≥0, ddyF(y)= α-24\,y^α/2-2 (1-y+y y )≥ 0, Moreover, F(1)=2α2(1−1)+α−22α⋅1⋅log⁡1−12⋅1⋅(1−1)=0.F(1)= 2α^2(1-1)+ α-22α· 1· 1- 12· 1·(1-1)=0. Hence, F(y)≥0F(y)≥ 0 for all y≥1y≥ 1, with strict inequality when y>1y>1 (equivalently, ϵ≠0ε≠ 0). When α<2α<2, we have ρ(ϵ;α,c)=2−α[(1+(ϵ/c)22−α)α/2−1].ρ(ε;α,c)= 2-α [ (1+ (ε/c)^22-α )^α/2-1 ]. Let y=1+(ϵ/c)22−α(≥1).y=1+ (ε/c)^22-α (≥ 1). Differentiating with respect to α gives ∂ρ∂α=−2α2(yα/2−1)+2−α2αyα/2log⁡y+12yα/2−1(y−1). ∂ρ∂α=- 2α^2 (y^α/2-1 )+ 2-α2α\,y^α/2 y+ 12\,y^α/2-1(y-1). For fixed α, ϵε, and c, define G:[1,∞)→ℝG:[1,∞) by G(y):=∂ρ∂α(ϵ,α,c)=−2α2(yα/2−1)+2−α2αyα/2log⁡y+12yα/2−1(y−1).G(y):= ∂ρ∂α(ε,α,c)=- 2α^2 (y^α/2-1 )+ 2-α2α\,y^α/2 y+ 12\,y^α/2-1(y-1). Differentiating in y gives dGdy=2−α4yα/2−2(ylog⁡y−(y−1)). dGdy= 2-α4\,y^α/2-2 (y y-(y-1) ). The fact that ylog⁡y−(y−1)≥0y y-(y-1)≥ 0 for y≥1y≥ 1 (a standard inequality, with equality only at y=1y=1) implies dGdy≥0 dGdy≥ 0, since 2−α>02-α>0 and yα/2−2≥0y^α/2-2≥ 0. As in the previous case, G(1)=0G(1)=0, so G(y)≥0G(y)≥ 0 for all y≥1y≥ 1, with strict inequality when y>1y>1 (i.e., ϵ≠0ε≠ 0). Combining both cases, we conclude that ∂ρ/∂α≥0∂ρ/∂α≥ 0 for all α≠0,2α≠ 0,2, with strict inequality unless ϵ=0ε=0. At α=0α=0 and α=2α=2, the function is not differentiable in α due to the absolute value, but the one-sided limits from the analysis above confirm that the loss increases as α moves away from these points whenever ϵ≠0ε≠ 0. ∎ F.2 Basic properties of the score function Lemma F.1. Fix α∈(0,1]α∈(0,1] and k=2−α∈[1,2)k=2-α∈[1,2). There exist constants Lα,Cα,c¯α>0L_α,C_α, c_α>0 depending only on α such that for all ε∈ℝ and c>0c>0: 1. (Bounded score) |ψ(ε;α,c)|≤min⁡|ε|c2,Lαc. |ψ( ;α,c) |≤ \ | |c^2,\ L_αc \. 2. (Derivative bound) Let ψ′(ε;α,c):=∂εψ(ε;α,c)ψ ( ;α,c):= ∂ ψ( ;α,c). Then |ψ′(ε;α,c)|≤Cαc2,ψ′(0;α,c)=1c2. |ψ ( ;α,c) |≤ C_αc^2, ψ (0;α,c)= 1c^2. 3. (Local positive curvature) If |ε|≤c | |≤ c, then ψ′(ε;α,c)≥c¯αc2.ψ ( ;α,c)≥ c_αc^2. Proof. Write u:=(ε/c)2/k≥0u:=( /c)^2/k≥ 0 and s:=1+u≥1s:=1+u≥ 1. Then ψ(ε;α,c)=εc2sα/2−1ψ( ;α,c)= c^2s^α/2-1. Since α/2−1≤0α/2-1≤ 0 and s≥1s≥ 1, we have sα/2−1≤1s^α/2-1≤ 1, hence |ψ|≤|ε|/c2 |ψ |≤ | |/c^2. Moreover, |ψ(ε;α,c)| |ψ( ;α,c) | decays for large |ε| | | (because α−1≤0α-1≤ 0), hence supε|ψ(ε;α,c)|≤Lα/c _ |ψ( ;α,c) |≤ L_α/c for some constant LαL_α. Differentiate to obtain ψ′(ε;α,c)=1c2sα/2−2[1+(α−1)u].ψ ( ;α,c)= 1c^2s^α/2-2 [1+(α-1)u ]. This yields ψ′(0)=1/c2ψ (0)=1/c^2 and |ψ′|≤Cα/c2 |ψ |≤ C_α/c^2. If |ε|≤c | |≤ c then u≤1/k≤1u≤ 1/k≤ 1 and 1+(α−1)u≥α1+(α-1)u≥α; also sα/2−2≥(1+1/k)α/2−2s^α/2-2≥(1+1/k)^α/2-2. Thus ψ′(ε;α,c)≥c¯α/c2ψ ( ;α,c)≥ c_α/c^2 with c¯α:=α(1+1/k)α/2−2>0 c_α:=α(1+1/k)^α/2-2>0. ∎ F.3 Proof of Theorem 5.3 We prove Theorem 5.3 by comparing the empirical score g^m g_m to its population counterpart gmg_m in a neighborhood of μ. The key idea is that, under finite variance and with cm→∞c_m→∞, the score behaves locally like the linear score of the sample mean, while still being well-defined even when α<1α<1. We first show that gmg_m identifies a unique population root θm _m near μ and that θm _m is close enough to μ so that its contribution is asymptotically negligible after m m scaling. We then show that g^m g_m has a unique central root with high probability and finally derive a Bahadur-type linearization around θm _m, which reduces the problem to a triangular-array CLT. Lemma F.2 (Population identification near μ). Assume Assumption 5.2 and fix α∈(0,1]α∈(0,1]. Let cm→∞c_m→∞. Then there exist constants c0>0c_0>0 and m0m_0 such that for all m≥m0m≥ m_0 and all x with |x−μ|≤cm/4 |x-μ |≤ c_m/4, −2cm2≤gm′(x)≤−c0cm2.- 2c_m^2\ ≤\ g_m (x)\ ≤\ - c_0c_m^2. In particular, gmg_m is strictly decreasing on x:|x−μ|≤cm/4\x: |x-μ |≤ c_m/4\, hence the equation gm(x)=0g_m(x)=0 has a unique solution on this interval, denoted by θm _m. Proof. Recall gm′(x)=−[ψ′(X−x;α,cm)]g_m (x)=-E[ψ (X-x;α,c_m)]. The upper bound follows from Lemma F.1(i). For the lower bound, define Ax:=|X−x|≤cmA_x:=\ |X-x |≤ c_m\. If |x−μ|≤cm/4 |x-μ |≤ c_m/4 then |X−μ|≤cm/2⊆Ax\ |X-μ |≤ c_m/2\ A_x, hence by Chebyshev, ℙ(Ax)≥1−4σ2cm2.P(A_x)≥ 1- 4σ^2c_m^2. On AxA_x, Lemma F.1(i) gives ψ′(X−x;α,cm)≥c¯α/cm2ψ (X-x;α,c_m)≥ c_α/c_m^2; on AxcA_x^c, Lemma F.1(i) gives |ψ′|≤Cα/cm2 |ψ |≤ C_α/c_m^2. Therefore [ψ′(X−x;α,cm)]≥c0/cm2E[ψ (X-x;α,c_m)]≥ c_0/c_m^2 for m large enough, implying gm′(x)≤−c0/cm2g_m (x)≤-c_0/c_m^2. ∎ Lemma F.3 (Population bias). Assume Assumption 5.2 and let θm _m be the unique solution of gm(θ)=0g_m(θ)=0 in |θ−μ|≤cm/4\ |θ-μ |≤ c_m/4\. Then |θm−μ|≤Cσ2cm, | _m-μ |≤ Cσ^2c_m, for a constant C depending only on α. Proof. The goal is to bound the (small) mismatch gm(μ)g_m(μ) induced by the robust score at finite cmc_m. Let ε:=X−μ :=X-μ and s:=1+(ε/cm)2/ks:=1+( /c_m)^2/k. Then ψ(ε;α,cm)=εcm2+εcm2(sα/2−1−1).ψ( ;α,c_m)= c_m^2+ c_m^2 (s^α/2-1-1 ). Since [ε]=0E[ ]=0, [ψ(X−μ;α,cm)]=1cm2[ε(sα/2−1−1)].E[ψ(X-μ;α,c_m)]= 1c_m^2E [ (s^α/2-1-1 ) ]. Split B=|ε|≤cmB=\ | |≤ c_m\ and BcB^c. On B, with u=(ε/cm)2/k∈[0,1]u=( /c_m)^2/k∈[0,1], a Taylor remainder bound gives |sα/2−1−1|≤Cu |s^α/2-1-1 |≤ Cu; hence |ε||sα/2−1−1|≤C|ε|u=C|ε|3cm2≤Cε2cm. | | |s^α/2-1-1 |≤ C | |u=C | |^3c_m^2≤ C ^2c_m. On BcB^c, |sα/2−1−1|≤1 |s^α/2-1-1 |≤ 1 and |ε||ε|>cm≤ε2/cm | |1\ | |>c_m\≤ ^2/c_m, so the same bound holds. Therefore |[ψ(X−μ;α,cm)]|≤Cσ2cm3,i.e.,|gm(μ)|≤Cσ2cm3. |E[ψ(X-μ;α,c_m)] |≤ Cσ^2c_m^3, .e., |g_m(μ) |≤ Cσ^2c_m^3. By the mean value theorem, for some θ¯ θ between μ and θm _m, 0=gm(θm)=gm(μ)+(θm−μ)gm′(θ¯).0=g_m( _m)=g_m(μ)+( _m-μ)g_m ( θ). Lemma F.2 ensures |gm′(θ¯)|≥c0/cm2 |g_m ( θ) |≥ c_0/c_m^2, hence |θm−μ|≤|gm(μ)||gm′(θ¯)|≤Cσ2/cm3c0/cm2=C′σ2cm. | _m-μ |≤ |g_m(μ) | |g_m ( θ) |≤ Cσ^2/c_m^3c_0/c_m^2= C σ^2c_m. ∎ Proof of Theorem 5.3. Step 1 (reduce to a CLT around the population root). Let θm _m be the unique solution of gm(θ)=0g_m(θ)=0 in |θ−μ|≤cm/4 |θ-μ |≤ c_m/4. Lemma F.3 yields |θm−μ|≤Cσ2/cm | _m-μ |≤ Cσ^2/c_m, hence m|θm−μ|≤Cσ2mcm=Cσ2m−γ→0. m\, | _m-μ |≤ Cσ^2 mc_m=Cσ^2m^-γ→ 0. Thus it suffices to prove m(X~m−θm)⇒N(0,σ2) m( X_m- _m) N(0,σ^2) and then apply Slutsky. (Informally, θm _m is the “population correction” induced by robustness, and it vanishes under our finite-variance tuning.) Step 2 (existence/uniqueness of the central root). We show that g^m g_m is strictly monotone on the central interval and that it changes sign at the endpoints. Using Lemma F.1(i) and Lemma F.2, one shows that with probability tending to one, infx∈[μ−rm,μ+rm]g^m′(x)≥c02cm2. _x∈[μ-r_m,μ+r_m] g_m (x)\ ≥\ c_02c_m^2. Hence g^m g_m is strictly monotone on [μ−rm,μ+rm][μ-r_m,μ+r_m]. Moreover, |ψ(⋅;α,cm)|≤Lα/cm |ψ(·;α,c_m) |≤ L_α/c_m implies supx∈[μ−rm,μ+rm]|g^m(x)−gm(x)|=Op(1cmm). _x∈[μ-r_m,μ+r_m] | g_m(x)-g_m(x) |=O_p\! ( 1c_m m ). Since gmg_m has slope ≍1/cm2 1/c_m^2 on this interval and rm/cm2≫1/(cmm)r_m/c_m^2 1/(c_m m) under our tuning, the endpoint signs of g^m g_m match those of gmg_m with probability →1→ 1, so g^m g_m crosses 0 exactly once on the interval. This yields the existence and uniqueness of X~m X_m on an event ℰmE_m with ℙ(ℰm)→1P(E_m)→ 1. Step 3 (Bahadur representation). On ℰmE_m, by the mean value theorem there exists X¯m X_m between X~m X_m and θm _m such that 0=g^m(X~m)=g^m(θm)+(X~m−θm)g^m′(X¯m),0= g_m( X_m)= g_m( _m)+( X_m- _m) g_m ( X_m), so X~m−θm=−g^m(θm)g^m′(X¯m). X_m- _m=- g_m( _m) g_m ( X_m). Since gm(θm)=0g_m( _m)=0, g^m(θm)=−1m∑i=1m[ψ(Xi−θm;α,cm)−ψ(X−θm;α,cm)]. g_m( _m)=- 1m _i=1^m [ψ(X_i- _m;α,c_m)-Eψ(X- _m;α,c_m) ]. Also g^m′(X¯m)=[ψ′(X−θm;α,cm)]+op(1/cm2)∼1/cm2 g_m ( X_m)=E[ψ (X- _m;α,c_m)]+o_p(1/c_m^2) 1/c_m^2, hence 1g^m′(X¯m)=cm2(1+op(1)). 1 g_m ( X_m)=c_m^2(1+o_p(1)). Therefore m(X~m−θm)=cm2m∑i=1m[ψ(Xi−θm;α,cm)−ψ(X−θm;α,cm)]+op(1). m( X_m- _m)= c_m^2 m _i=1^m [ψ(X_i- _m;α,c_m)-Eψ(X- _m;α,c_m) ]+o_p(1). This linearization makes explicit that the leading term is an average of i.i.d. centered scores. Step 4 (triangular-array CLT). Let Zi,m:=cm2[ψ(Xi−θm;α,cm)−ψ(X−θm;α,cm)]Z_i,m:=c_m^2 [ψ(X_i- _m;α,c_m)-Eψ(X- _m;α,c_m) ]. Then cm2ψ(X−θm;α,cm)→X−μc_m^2ψ(X- _m;α,c_m)→ X-μ a.s., and |cm2ψ(⋅)|≤|⋅| |c_m^2ψ(·) |≤ |· |; dominated convergence yields Var(Zi,m)→σ2Var(Z_i,m)→σ^2. Lindeberg holds because [(X−μ)2]<∞E[(X-μ)^2]<∞. Thus 1m∑i=1mZi,m⇒N(0,σ2) 1 m _i=1^mZ_i,m N(0,σ^2). Combining with Step 3 gives m(X~m−θm)⇒N(0,σ2) m( X_m- _m) N(0,σ^2), and Step 1 finishes. ∎ F.4 Proof of Theorem 5.5 The proof follows the same three-step structure as the finite-variance case: (i) identify a unique population root θm _m and bound its bias relative to μ; (i) control the empirical score at θm _m with high probability; and (i) convert score control into an error bound for the empirical root via a mean-value expansion. The only difference is that, under a mere (1+ϵ)(1+ε)-moment, the scale cmc_m must be tuned at the deviation level to balance truncation bias and concentration. Lemma F.4 (Bias under (1+ϵ)(1+ε)-moment). Assume Assumption 5.4 and fix α∈(0,1]α∈(0,1]. Let θm _m be the unique root of gm(θ)=0g_m(θ)=0 in |θ−μ|≤cm/4\ |θ-μ |≤ c_m/4\. Then |θm−μ|≤Cv1+ϵcmϵ. | _m-μ |≤ C v_1+εc_m^ε. Proof. The argument parallels Lemma F.3, replacing second-moment tail bounds by (1+ϵ)(1+ε)-moment bounds (Markov). One shows |gm(μ)|≤Cv1+ϵ/cm2+ϵ |g_m(μ) |≤ Cv_1+ε/c_m^2+ε and |gm′(⋅)|≍1/cm2 |g_m (·) | 1/c_m^2 on |x−μ|≤cm/4 |x-μ |≤ c_m/4, yielding |θm−μ|≤Cv1+ϵ/cmϵ | _m-μ |≤ Cv_1+ε/c_m^ε. ∎ Lemma F.5 (Variance bound for ψ). Assume Assumption 5.4. Then for any x∈ℝx and c>0c>0, Var(ψ(X−x;α,c))≤Cv1+ϵc3+ϵ,Var (ψ(X-x;α,c) )≤ C v_1+εc^3+ε, for a constant C depending only on α and ϵε. Proof. By Lemma F.1(i), ψ(ε;α,c)2≤ε2c4|ε|≤c+Lα2c2|ε|>c.ψ( ;α,c)^2≤ ^2c^41\ | |≤ c\+ L_α^2c^21\ | |>c\. For the first term, use |ε|2|ε|≤c≤c1−ϵ|ε|1+ϵ | |^21\ | |≤ c\≤ c^1-ε | |^1+ε. For the second term, Markov gives ℙ(|ε|>c)≤[|ε|1+ϵ]/c1+ϵP( | |>c) [ | |^1+ε]/c^1+ε. Combine these bounds and use [|X−x|1+ϵ]≤C([|X−μ|1+ϵ]+|x−μ|1+ϵ)E[ |X-x |^1+ε]≤ C(E[ |X-μ |^1+ε]+ |x-μ |^1+ε). ∎ Lemma F.6 (Bernstein control). Assume Assumption 5.4 and let θm _m be the population root. Then for any δ∈(0,1/2)δ∈(0,1/2), with probability at least 1−δ1-δ, |g^m(θm)|≤C(v1+ϵlog⁡(2/δ)mcm3+ϵ+log⁡(2/δ)m⋅1cm). | g_m( _m) |≤ C ( v_1+ε (2/δ)m\,c_m^3+ε+ (2/δ)m· 1c_m ). Proof. We control the centered empirical score by a bounded-variance Bernstein inequality. Let Yi:=ψ(Xi−θm;α,cm)−[ψ(X−θm;α,cm)]Y_i:=ψ(X_i- _m;α,c_m)-E[ψ(X- _m;α,c_m)]. Then g^m(θm)=−(1/m)∑i=1mYi g_m( _m)=-(1/m) _i=1^mY_i. Lemma F.1(i) gives |Yi|≤2Lα/cm |Y_i |≤ 2L_α/c_m. Lemma F.5 gives Var(Yi)≤Cv1+ϵ/cm3+ϵVar(Y_i)≤ Cv_1+ε/c_m^3+ε. Apply Bernstein’s inequality for bounded independent variables. ∎ Proof of Theorem 5.5. Let θm _m be the population root. We first relate the estimation error to the empirical score at θm _m. By the mean value theorem (on the event where g^m g_m is strictly monotone and crosses zero in ℐmI_m), X~m−θm=−g^m(θm)g^m′(X¯m) X_m- _m=- g_m( _m) g_m ( X_m) for some X¯m X_m between X~m X_m and θm _m. The identification argument (population [ψ′]≍1/cm2E[ψ ] 1/c_m^2 and a uniform LLN for g^m′ g_m on ℐmI_m) gives g^m′(X¯m)≥c02cm2 g_m ( X_m)≥ c_02c_m^2 with probability at least 1−δ/21-δ/2 for m large enough. Hence |X~m−θm|≤Ccm2|g^m(θm)|. | X_m- _m |≤ Cc_m^2 | g_m( _m) |. Next, apply Lemma F.6 (with probability at least 1−δ/21-δ/2) to obtain |X~m−θm|≤C(v1+ϵlog⁡(2/δ)mcm1−ϵ2+log⁡(2/δ)mcm). | X_m- _m |≤ C ( v_1+ε (2/δ)m\,c_m 1-ε2+ (2/δ)m\,c_m ). Finally, we add the population bias controlled by Lemma F.4: |X~m−μ|≤|X~m−θm|+|θm−μ|≤C(v1+ϵlog⁡(2/δ)mcm1−ϵ2+log⁡(2/δ)mcm+v1+ϵcmϵ). | X_m-μ |≤ | X_m- _m |+ | _m-μ |≤ C ( v_1+ε (2/δ)m\,c_m 1-ε2+ (2/δ)m\,c_m+ v_1+εc_m^ε ). The chosen tuning cm=cm(δ)=τ(v1+ϵmlog⁡(2/δ))11+ϵc_m=c_m(δ)=τ ( v_1+εm (2/δ) ) 11+ε balances the bias term v1+ϵ/cmϵv_1+ε/c_m^ε with the leading stochastic term, and a direct calculation shows that each term on the right-hand side is bounded by C1(v1+ϵlog⁡(C2/δ)m)ϵ1+ϵC_1 ( v_1+ε (C_2/δ)m ) ε1+ε for suitable constants C1,C2C_1,C_2. ∎ Appendix G Proofs for the ARE results G.1 A median-boosting lemma Lemma G.1 (Median boosting). Let Z1,…,ZkZ_1,…,Z_k be independent random variables and let θ∈ℝθ . Assume that for some r>0r>0 and p∈(0,1/2)p∈(0,1/2), ℙ(|Zj−θ|≤r)≥1−p,j=1,…,k.P (|Z_j-θ|≤ r )≥ 1-p, j=1,…,k. Let median⁡(Z1,…,Zk)median(Z_1,…,Z_k) be any median. Then ℙ(|median⁡(Z1,…,Zk)−θ|>r)≤exp⁡(−2k(1/2−p)2).P\! ( |median(Z_1,…,Z_k)-θ |>r )≤ \! (-2k(1/2-p)^2 ). In particular, if p=1/4p=1/4, then ℙ(|median⁡(Z1,…,Zk)−θ|>r)≤e−k/8P(|median(Z_1,…,Z_k)-θ|>r)≤ e^-k/8. Proof. Let Ij:=|Zj−θ|≤rI_j:=1\|Z_j-θ|≤ r\. Then [Ij]≥1−p>1/2E[I_j]≥ 1-p>1/2 and Ij\I_j\ are independent. If |median⁡(Z1,…,Zk)−θ|>r|median(Z_1,…,Z_k)-θ|>r, then at least ⌈k/2⌉ k/2 of the ZjZ_j’s must fall outside [θ−r,θ+r][θ-r,θ+r], i.e., ∑j=1kIj≤k/2 _j=1^kI_j≤ k/2. Hence, by Hoeffding’s inequality, ℙ(∑j=1kIj≤k/2)=ℙ(∑j=1k(Ij−Ij)≤−k(I1−1/2))≤exp⁡(−2k(I1−1/2)2)≤exp⁡(−2k(1/2−p)2).P\! ( _j=1^kI_j≤ k/2 )=P\! ( _j=1^k(I_j-EI_j)≤-k(EI_1-1/2) )≤ \! (-2k(EI_1-1/2)^2 )≤ \! (-2k(1/2-p)^2 ). ∎ G.2 Proof of Theorem 5.7 Proof of Theorem 5.7. Fix δ0:=1/4 _0:=1/4 and let m=⌊n/k⌋m= n/k . Apply Theorem 5.5 on each block (with confidence level δ0 _0) to obtain: there exist constants C1′,C2′>0C_1 ,C_2 >0 such that for all sufficiently large m, ℙ(|X~j−μ|≤rm:=C1′(v1+ϵlog⁡(C2′/δ0)m)ϵ1+ϵ)≥1−δ0=34.P\! (| X_j-μ|≤ r_m:=C_1 ( v_1+ε (C_2 / _0)m ) ε1+ε )≥ 1- _0= 34. Since the blocks are disjoint, X~1,…,X~k X_1,…, X_k are independent. Therefore Lemma G.1 (with p=δ0=1/4p= _0=1/4) yields ℙ(|μ~−μ|>rm)≤e−k/8.P (| μ-μ|>r_m )≤ e^-k/8. Choosing k=⌈8log⁡(2/δ)⌉k= 8 (2/δ) gives e−k/8≤δ/2e^-k/8≤δ/2. It remains to rewrite rmr_m in terms of n and δ. Since m=⌊n/k⌋≥n/(2k)m= n/k ≥ n/(2k) for n large enough, rm≤C(v1+ϵkn)ϵ1+ϵ≤C(v1+ϵlog⁡(2/δ)n)ϵ1+ϵ.r_m≤ C ( v_1+ε\,kn ) ε1+ε≤ C ( v_1+ε (2/δ)n ) ε1+ε. Absorb constants and the log⁡(2/δ) (2/δ) vs log⁡(C2/δ) (C_2/δ) adjustment into C1,C2C_1,C_2 to conclude the claim. ∎ G.3 Proof of Theorem 5.6 Lemma G.2 (Median CLT with local density). Let Y1,…,YkY_1,…,Y_k be i.i.d. with distribution function F and density f. Assume F(μ)=1/2F(μ)=1/2, f(μ)>0f(μ)>0, and f is continuous at μ. Let μ^:=median⁡(Y1,…,Yk) μ:=median(Y_1,…,Y_k). Then k(μ^−μ)⇒N(0,14f(μ)2). k\,( μ-μ) N\! (0, 14f(μ)^2 ). Proof. This is the classical sample-median CLT (e.g., via the Bahadur representation for sample quantiles). ∎ Proof of Theorem 5.6. Let k=kn→∞k=k_n→∞ and m=⌊n/k⌋→∞m= n/k →∞. By Theorem 5.3, for each block, m(X~j−μ)⇒N(0,σ2),j=1,…,k. m\,( X_j-μ) N(0,σ^2), j=1,…,k. Let FmF_m and fmf_m denote the cdf and density of X~j X_j. Under the additional local normal approximation condition stated in the theorem, we have Fm(μ)=1/2F_m(μ)=1/2 and fm(μ)=mσ2π(1+o(1)).f_m(μ)= mσ 2π\,(1+o(1)). Apply Lemma G.2 to Yj=X~jY_j= X_j with F=FmF=F_m and f=fmf=f_m: k(μ~−μ)⇒N(0,14fm(μ)2)=N(0,π2⋅σ2m). k\,( μ-μ) N\! (0, 14f_m(μ)^2 )=N\! (0, π2· σ^2m ). Multiplying both sides by m m gives mk(μ~−μ)=n(μ~−μ)⇒N(0,π2σ2). mk\,( μ-μ)= n\,( μ-μ) N\! (0, π2σ^2 ). Finally, since X¯n X_n has asymptotic variance σ2/nσ^2/n, we obtain ARE(μ~,X¯n)=2/πARE( μ, X_n)=2/π. ∎ Appendix H Hyperparameter Settings for The Experiments Table 9: Training hyperparameters for dual-agent ILR training. Hyperparameter Value Description Learning rate (Agent 1) 1×10−61× 10^-6 Adam optimizer Learning rate (Agent 2) 1×10−61× 10^-6 Adam optimizer Sampling temperature 0.50.5 — Training batch size 3232 — Number of episodes 1010 — Max prompt length 61446144 Tokens Max generation length 15361536 Tokens Top-p sampling 0.950.95 Nucleus sampling Discount factor γ 1.01.0 — KL estimator K2 — Advantage estimator ILR — Precision BF16 Mixed precision Table 10: Reward noise hyperparameters (Sparse Cauchy Outlier Noise). Parameter Symbol Value Description Noise type — Outlier Sparse one-sided heavy-tail Group contamination prob. p 0.20.2 Probability of corrupting a group Cauchy scale γ 1.01.0 Controls distribution width Spike scale α 10.010.0 Noise magnitude scaling factor Spike clip Δmax _ 10.010.0 Maximum absolute noise value Noise target — Argmax Inject noise to group maximum All experiments are conducted on a single node with 8 A800 GPUs. We use Ray for distributed training orchestration and vLLM for efficient inference with 4 inference engines and tensor parallelism of 2. Gradient checkpointing and Adam parameter offloading are enabled to optimize GPU memory usage. To simulate reward uncertainty in realistic scenarios, we introduce a sparse one-sided heavy-tail spike noise based on the Cauchy distribution.

Instant research discovery