Paper deep dive
HO-FL: Hybrid-Order Federated Learning for Heterogeneous Edge Devices
Qiyuan Chen, Xian Wu, Yanan Ma, Xianhao Chen
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Federated learning (FL) on memory-constrained edge devices faces a dilemma: first-order (FO) optimization (i.e., backpropagation) demands substantial memory, whereas zeroth-order (ZO) optimization suffers from severe convergence slowdown. To resolve this dilemma, we introduce HO-FL, a hybrid-order FL framework that trains a model's bottom segment with ZO optimization and its top segment with FO optimization. Each device can flexibly select its order boundary according to its memory budget while participating in the training of the same global model. Moreover, our convergence analysis reveals a new, fundamental trade-off: clients with larger FO-trained segments can provide more accurate updates, but favoring them can underrepresent other clients' data. We connect this trade-off to the bias and variance of actual multi-step local updates, yielding a sampling optimization problem and a practical dimension-aware approximation with direct model averaging. Experiments on language tasks examine task performance, client memory, and sampling under data heterogeneity. The results show that hybrid-order local training can retain much of the full-FO performance with substantially lower client memory requirements. Our code is available at this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2609.39074v1
- Canonical: https://arxiv.org/abs/2609.39074v1
Trouble viewing inline? Open PDF directly →
Full Text
109,952 characters extracted from source content.
Expand or collapse full text
HO-FL: Hybrid-Order Federated Learning for Heterogeneous Edge Devices Qiyuan Chen †thanks: Equal contribution. Affiliation: University of Hong Kong (hku.hk) Email: qiyuanchen@connect.hku.hk Xian Wu Affiliation: University of Hong Kong (hku.hk) Email: u3012772@connect.hku.hk Yanan Ma Affiliation: City University of Hong Kong (cityu.edu.hk) Email: yananma8-c@my.cityu.edu.hk Xianhao Chen †thanks: Corresponding author: Xianhao Chen (xcheneee@hku.hk). Affiliation: University of Hong Kong (hku.hk) Email: xcheneee@hku.hk Abstract Federated learning (FL) on memory-constrained edge devices faces a dilemma: first-order (FO) optimization (i.e., backpropagation) demands substantial memory, whereas zeroth-order (ZO) optimization suffers from severe convergence slowdown. To resolve this dilemma, we introduce HO-FL, a hybrid-order FL framework that trains a model’s bottom segment with ZO optimization and its top segment with FO optimization. Each device can flexibly select its order boundary according to its memory budget while participating in the training of the same global model. Moreover, our convergence analysis reveals a new, fundamental trade-off: clients with larger FO-trained segments can provide more accurate updates, but favoring them can underrepresent other clients’ data. We connect this trade-off to the bias and variance of actual multi-step local updates, yielding a sampling optimization problem and a practical dimension-aware approximation with direct model averaging. Experiments on language tasks examine task performance, client memory, and sampling under data heterogeneity. The results show that hybrid-order local training can retain much of the full-FO performance with substantially lower client memory requirements. Our code is available at https://github.com/HKU-WILL-Lab/HO-FL. 1 Introduction Deploying and fine-tuning large language models (LLMs) directly on edge devices provides a promising paradigm for unlocking decentralized real-world data while preserving user data privacy (Wu et al., 2026). Federated learning (FL) (McMahan et al., 2017) has been widely explored as an effective framework to collaboratively adapt models across distributed clients without centralizing raw data. However, standard FL relies on first-order (FO) optimization via backpropagation (BP), which requires storing intermediate activations across model layers. For modern LLMs, these activations easily exceed the physical memory limits of commodity edge hardware, posing a severe memory bottleneck for on-device federated adaptation. To mitigate this activation memory overhead, zeroth-order (ZO) optimization has been increasingly adopted in federated edge learning (Fang et al., 2022; Li et al., 2025). By estimating gradient directions solely through perturbed forward passes, ZO optimization avoids constructing a backward computation graph, thereby reducing client training memory to the inference level. Nevertheless, pure ZO gradient optimization incurs a severe variance penalty: random directional perturbations produce gradient estimators whose variance scales proportionally with the trainable parameter dimension d. Consequently, pure ZO federated algorithms typically require orders-of-magnitude more local iterations and communication rounds than their FO counterparts to reach target performance (Chen et al., 2026; Li et al., 2026), limiting their practicality for complex downstream tasks. Recent efforts on hybrid-order (HO) optimization provide an appealing methodology to overcome this memory–accuracy dilemma. By partitioning the network architecture, HO optimization structurally decouples the optimization landscape: the bottom layers are optimized via BP-free ZO perturbations, while the top layers are trained using FO gradients. While early exploratory efforts have leveraged such hybrid updates within split learning paradigms (Chen et al., 2026; Lnu et al., 2026; Kou et al., 2026), directly adapting hybrid-order optimization to practical FL is hindered by heterogeneous edge constraints. In federated settings, real-world edge devices exhibit diverse physical memory capacities; enforcing a globally identical order boundary either exceeds the memory budget of resource-constrained devices or underutilizes the capacity of more capable ones. These challenges lead to a natural question: Can we design a fully federated hybrid-order optimization paradigm that accommodates heterogeneous edge hardware while retaining the convergence advantages of first-order training? To answer this question, we propose HO-FL, a hybrid-order FL framework tailored for heterogeneous edge devices. In HO-FL, each client executes both ZO updates on the bottom segment and FO backpropagation on the top segment entirely on-device, performing multi-step local training and synchronizing with the central server only at the end of each round. To accommodate hardware disparity, HO-FL enables each client to configure a personalized order boundary matched to its local memory budget. To alleviate uplink communication overhead, HO-FL further encodes bottom ZO updates into compact scalar-seed pairs, requiring full-vector uploads only for the top FO parameters. We establish the non-convex convergence of HO-FL, proving that the average gradient norm is bounded by ((1+1q∑i=1NpidiZO)/T)+(ζ2N∑i=1N(pi−1N)2)O ( (1+ 1q _i=1^Np_id_i^ZO)/T )+O (ζ^2N _i=1^N(p_i- 1N)^2 ). Here, the first term characterizes the optimization rate governed by the sampling-weighted ZO parameter dimension, while the second term quantifies the representation bias induced by client data heterogeneity ζ2ζ^2. Furthermore, our bound provides a unified theoretical guarantee, seamlessly recovering existing hybrid-order and federated convergence rates as special cases when clients participate uniformly (pi=1/Np_i=1/N) or share homogeneous order boundaries. Crucially, this bound reveals a promising optimization opportunity: the two terms formalize a fundamental trade-off between local update quality and global data representation. While uniform sampling is unavoidably bottlenecked by noisy updates from ZO-dominated clients, non-uniform selection p can actively suppress gradient variance by prioritizing FO-dominant clients at the cost of a controlled representation penalty. Grounded in this trade-off, we design an efficient Dimension-Aware Client Sampling policy via dependent rounding, provably and empirically outperforms uniform sampling. Our main contributions are summarized as follows: • Framework & Unified Convergence Theory: We propose HO-FL, the first fully federated hybrid-order optimization framework supporting client-adaptive order boundaries and compact seed-based uplink transmission. We establish a unified convergence theory for HO-FL, showing that the convergence bound decomposes into an optimization term ((1+1q∑i=1NpidiZO)/T)O ( (1+ 1q _i=1^Np_id_i^ZO)/T ) and a sampling bias term (ζ2N∑i=1N(pi−1N)2)O (ζ^2N _i=1^N(p_i- 1N)^2 ), seamlessly subsuming existing hybrid-order analyses as special cases. • Dimension-Aware Client Sampling: Grounded in our convergence theory, we formulate a dimension-aware client sampling policy via dependent rounding. Our strategy selectively prioritizes low-variance updates while bounding sampling bias, which is demonstrated both theoretically and empirically to accelerate global convergence over uniform client sampling. • Empirical Validation: We evaluate HO-FL across diverse language understanding and generation benchmarks on OPT-125M, Qwen2.5-1.5B, and SmolLM3-3B. The results demonstrate that HO-FL closely approaches full first-order federated performance while slashing client peak memory by up to 77.4% on the most resource-constrained devices and substantially outperforming pure ZO baselines. 2 Related Work Zeroth-order optimization in federated learning. Zeroth-order (ZO) optimization avoids activation caching by estimating gradients via randomized finite differences, reducing client memory to inference levels (Malladi et al., 2023; Fang et al., 2022). To curb the communication overhead of exchanging perturbations, DeComFL (Li et al., 2025) and HiSo (Li et al., 2026) employ pseudo-random seed synchronization and Hessian guidance, respectively. Nevertheless, pure ZO federated algorithms inevitably suffer from gradient estimation variance that scales with the full parameter dimension d, severely slowing convergence. HO-FL circumvents this bottleneck by confining ZO optimization strictly to memory-critical bottom segments. Hybrid-order optimization across distributed paradigms. To balance memory footprint and gradient accuracy, recent works partition networks into a bottom ZO segment and a top FO segment. Although existing methods realize this hybrid structure differently, they can be cast into a unified framework governed by boundary activations and localized surrogate objectives (see Appendix A for details): HOSL (Lnu et al., 2026) evaluates perturbations through end-to-end forward passes; HO-SFL (Chen et al., 2026) freezes cut-layer activation feedback to construct a linear surrogate; and HERON-SFL (Kou et al., 2026) introduces a client-side auxiliary network to decouple perturbations from server dependencies. However, these methods are fundamentally tied to split learning (SL) architectures, which require per-iteration activation transmissions and enforce identical order boundaries across all clients. In contrast, HO-FL executes both segments fully locally, allows client-adaptive order boundaries to fit their heterogeneous memory limits, and enables dimension-aware client sampling to enhance performance. 3 HO-FL: Hybrid-Order Federated Optimization Problem formulation. We consider an FL system with N distributed edge clients, where each client i∈1,…,Ni∈\1,…,N\ holds a local data distribution iD_i. Given global model parameters w∈ℝdw ^d and a sample loss ℓ(w,x,y) (w;x,y), the global objective is: minw∈ℝdf(w)≔1N∑i=1Nfi(w),wherefi(w)≔(x,y)∼i[ℓ(w,x,y)]. _w ^df(w) 1N _i=1^Nf_i(w), f_i(w) _(x,y) _i [ (w;x,y) ]. (1) Client sampling via dependent rounding. At each global round t, the server samples a cohort StS_t of K distinct clients among N clients to participate in training. We characterize client participation via normalized selection probabilities pip_i, which reflect the long-term proportion of selections allocated to client i over the cumulative sampling budget, satisfying pi∈[0,1/K]p_i∈[0,1/K] and ∑i=1Npi=1 _i=1^Np_i=1. While multiple sampling schemes can realize this marginal distribution, we draw the K distinct participants via dependent rounding (Gandhi et al., 2006), implemented as St=DepRound(Kp)S_t=DepRound(Kp) in Algorithm 2. This choice ensures theoretical tractability, with structural properties detailed in Appendix C.4. The server then dispatches independent pseudorandom seeds seedi,t,jee=0,j=1E−1,q\seed_i,t,j^e\_e=0,j=1^E-1,q to each selected client i∈Sti∈ S_t for subsequent local ZO updates. (a) (b) Figure 1: Overview of the HO-FL framework. (a) Client-side hybrid-order step: the bottom segment fi,ZOf_i,ZO is updated via ZO perturbations evaluated on surrogate objective ℒiL_i guided by order-boundary activation feedback δi _i, while the top segment fi,FOf_i,FO performs standard BP. (b) System-level federated workflow: the server selects cohorts via dimension-aware dependent rounding, broadcasts wtw_t, and reconstructs full model updates in ℝdR^d via seed replay for global model aggregation. Memory-adaptive boundaries. Client i chooses a boundary bib_i according to its memory constraint and partitioning its trainable parameters into two parts: wi=[(wiZO)⊤,(wiFO)⊤]⊤,w_i= [(w_i^ZO) ,(w_i^FO) ] , (2) where wiZO∈ℝdiZOw_i^ZO ^d_i^ZO and wiFO∈ℝdiFOw_i^FO ^d_i^FO denote the bottom and top parameter segments, respectively, satisfying diZO+diFO=d_i^ZO+d_i^FO=d. Since clients can update the ZO segment without caching the corresponding activation, flexible order boundaries allow resource-constrained devices to participate by expanding wiZOw_i^ZO, while enabling more capable devices to retain larger FO segments. A hybrid-order local step. As illustrated in Figure 1, we represent the neural network on client i via a bottom sub-network fi,ZO(⋅,wiZO)f_i,ZO(·;w_i^ZO) and a top sub-network fi,FO(⋅,wiFO)f_i,FO(·;w_i^FO). For a minibatch ξ=(x,y)ξ=(x,y), the bottom segment first computes the boundary activation zi=fi,ZO(x,wiZO)z_i=f_i,ZO(x;w_i^ZO) without caching intermediate activations. Feeding ziz_i into the top segment yields the task prediction y^i=fi,FO(zi,wiFO) y_i=f_i,FO(z_i;w_i^FO) and its associated loss ℓ(y^i,y) ( y_i,y). Subsequently, backpropagation through the top segment produces both the parameter gradient with respect to wiFOw_i^FO and the activation gradient with respect to ziz_i: giFO=∇wiFOℓ(wi,ξ),δi=∇ziℓ(fi,FO(zi,wiFO),y).g_i^FO= _w_i^FO (w_i;ξ), _i= _z_i (f_i,FO(z_i;w_i^FO),y). (3) Following HO-SFL (Chen et al., 2026), we construct a local surrogate objective for the bottom parameters v: ℒi(v)=⟨δi,fi,ZO(x,v)⟩.L_i(v)= _i,f_i,ZO(x;v) . (4) By the chain rule, the gradient of this surrogate objective with respect to v exactly recovers the first-order gradient of the original loss. Introducing this surrogate eliminates the need for end-to-end forward passes with perturbed parameters, restricting ZO evaluations strictly to the bottom segment. Using pseudorandom seeds dispatched by the server, the client generates q independent random perturbation vectors via a random number generator (RNG): ui,j=RNG(seedi,j)∼(0,IdiZO),j=1,…,q.u_i,j=RNG(seed_i,j) (0,I_d_i^ZO), j=1,…,q. (5) The client then estimates the ZO gradient for the bottom segment using simultaneous perturbation stochastic approximation (SPSA): si,j=ℒi(wiZO+μui,j)−ℒi(wiZO)μ,g^iZO=1q∑j=1qsi,jui,j,s_i,j= L_i(w_i^ZO+μ u_i,j)-L_i(w_i^ZO)μ, g_i^ZO= 1q _j=1^qs_i,ju_i,j, (6) In our framework, each client performs this hybrid-order update for E local steps per communication round. Algorithm 3 combines the forward pass, top-segment BP, and bottom-segment estimation in a single local loop. The boundary gradient δi _i is recomputed at each step and held fixed within its q finite differences. Indexing the global round by t and the local iteration by e, both segment gradients are evaluated at wi,tew_i,t^e before the updates: wi,tZO,e+1=wi,tZO,e−ηg^i,tZO,e,wi,tFO,e+1=wi,tFO,e−ηgi,tFO,e,e=0,…,E−1. aligned w_i,t^ZO,e+1&=w_i,t^ZO,e-η g_i,t^ZO,e,\\ w_i,t^FO,e+1&=w_i,t^FO,e-η g_i,t^FO,e, aligned e=0,…,E-1. (7) Client uploads. Upon completing E local steps, client i uploads only the sequence of directional projection scalars si,t,je\s_i,t,j^e\ and the accumulated top-segment parameter difference: Δi,tFO=wi,tFO,E−wi,tFO,0. _i,t^FO=w_i,t^FO,E-w_i,t^FO,0. (8) Because the random seeds are pre-assigned by the server, the client does not need to transmit the high-dimensional perturbation vectors. Compared to conventional FL where full model deltas in ℝdR^d must be uploaded, the ZO component in HO-FL achieves dimension-free uplink transmission. As a consequence, configuring a larger ZO segment not only alleviates peak activation memory on the client, but also substantially reduces uplink communication overhead. Server reconstruction and aggregation. Upon collecting the local updates from the selected participant set StS_t, the server calls Reconstruct (Algorithm 4) to recover the bottom-segment parameter updates using the known seeds: ui,t,je=RNG(seedi,t,je),Δi,tZO=−ηq∑e=0E−1∑j=1qsi,t,jeui,t,je.u_i,t,j^e=RNG(seed_i,t,j^e), _i,t^ZO=- ηq _e=0^E-1 _j=1^qs_i,t,j^eu_i,t,j^e. (9) The server then concatenates Δi,tZO _i,t^ZO and the received Δi,tFO _i,t^FO according to client i’s order boundary bib_i, reconstructing the full model update Δi,t=[(Δi,tZO)⊤,(Δi,tFO)⊤]⊤∈ℝd _i,t=[( _i,t^ZO) ,( _i,t^FO) ] ^d. Although participating devices adopt heterogeneous order boundaries, their updates are fully compatible in the shared parameter space, enabling direct coordinate-wise model averaging: wt+1=wt+1K∑i∈StΔi,t.w_t+1=w_t+ 1K _i∈ S_t _i,t. (10) The updated model wt+1w_t+1 is subsequently broadcast to clients for the next round. Under this design, the client uplink transmission per round requires only (Eq+diFO)O(Eq+d_i^FO) elements, which scales gracefully even for overparameterized architectures. For stateful optimizers such as AdamW, the server exactly mirrors the client trajectory via deterministic replay, as elaborated in Appendix B.3. Algorithm 1 HO-FL: hybrid-order federated optimization 1: w0w_0; bi,ii=1N\b_i,D_i\_i=1^N; T,K,E,q,μ,η,βT,K,E,q,μ,η,β. 2: Set r by (16). 3: pi←(1−β)/N+βrip_i←(1-β)/N+β r_i, i=1,…,Ni=1,…,N. 4: for t=0,…,T−1t=0,…,T-1 do 5: // Server: sampling and broadcast 6: St←S_t← DepRound(Kp)(Kp). ⊳ Algorithm 2 7: Broadcast wtw_t; dispatch independent seedi,t,jee,j\seed_i,t,j^e\_e,j to each i∈Sti∈ S_t. 8: // Clients: local hybrid-order training 9: for each i∈Sti∈ S_t in parallel do 10: (si,t,jee,j,Δi,tFO)←(\s_i,t,j^e\_e,j, _i,t^FO)← ClientUpdate(i,t,wt)(i,t,w_t). ⊳ Algorithm 3 11: Upload si,t,jee,j\s_i,t,j^e\_e,j and Δi,tFO _i,t^FO. 12: end for 13: // Server: reconstruction and direct averaging 14: for each i∈Sti∈ S_t do 15: Δi,t← _i,t← Reconstruct(i,t,si,t,jee,j,Δi,tFO)(i,t,\s_i,t,j^e\_e,j, _i,t^FO). ⊳ Algorithm 4 16: end for 17: wt+1←wt+1K∑i∈StΔi,tw_t+1← w_t+ 1K _i∈ S_t _i,t. ⊳ Eq. (10) 18: end for 19: return wTw_T. 4 Convergence Analysis In this section, we analyze the theoretical convergence of the HO-FL local updates in (7) under non-convex objectives, establishing the theoretical foundation for our sampling policy in Section 5. 4.1 Theoretical Assumptions Our analysis builds upon standard regularity conditions in non-convex distributed and zeroth-order optimization (Karimireddy et al., 2020; Chen et al., 2026). Specifically, we assume that each client objective fif_i is L-smooth and the global loss is bounded below (Assumption C.1); reference stochastic gradients are conditionally unbiased with bounded variance (Assumption C.2); the local surrogate ℒiL_i is LisL_i^s-smooth (Assumption C.3); and inter-client data divergence satisfies the affine heterogeneity condition with residual heterogeneity bound ζ2ζ^2 (Assumption C.4). Formal mathematical statements and regularity discussions are deferred to Appendix C.1. 4.2 Non-Convex Convergence Rate By controlling the multi-step local trajectory drift and hybrid gradient estimation variance (with intermediate descent guarantees deferred to Appendix C.5), and properly calibrating the step size and perturbation radius, we establish the following convergence rate of HO-FL. Theorem 4.1 (Informal Statement of Non-Convex Convergence). Under Assumptions C.1–C.4, with the local step size η and perturbation radius μ suitably calibrated, the iterates of HO-FL satisfy: 1T∑t=0T−1‖∇f(wt)‖2≤(1+1q∑i=1NpidiZOT)⏟Optimization Rate+(ζ2N∑i=1N(pi−1N)2)⏟Sampling Representation Bias. 1T _t=0^T-1E\|∇ f(w_t)\|^2≤ O\! ( 1+ 1q _i=1^Np_id_i^ZOT )_Optimization Rate+ O\! (ζ^2N _i=1^N (p_i- 1N )^2 )_Sampling Representation Bias. (11) Here, ζ2ζ^2 represents the client data heterogeneity parameter defined in Assumption C.4. The formal statement with explicit step-size conditions and the complete proof is provided in Appendix C.6. Theoretical generality. Theorem 4.1 provides a unified theoretical guarantee that subsumes existing hybrid-order paradigms as special cases. Specifically, under uniform client participation (pi=1/Np_i=1/N) and a global identical order boundary (diZO≡dZOd_i^ZO≡ d^ZO), the sampling representation bias vanishes, reducing the bound in (11) to: 1T∑t=0T−1‖∇f(wt)‖2≤(1+dZO/qT). 1T _t=0^T-1E\|∇ f(w_t)\|^2 \! ( 1+d^ZO/qT ). (12) This rate exactly recovers the non-convex convergence established in prior hybrid-order frameworks such as HO-SFL (Chen et al., 2026). Furthermore, at the architectural extremes, Eq. (12) naturally collapses to the canonical (1/T)O(1/ T) rate of FedAvg (Wang et al., 2024) for full first-order training (dZO=0d^ZO=0), and the ((1+d/q)/T)O( (1+d/q)/T) rate of pure ZO-FL (Fang et al., 2022; Li et al., 2025) when updates are entirely zeroth-order (dZO=d^ZO=d). This demonstrates that our convergence bound offers a strictly more general theoretical framework. 5 Dimension-Aware Client Sampling Theorem 4.1 exposes an optimization opportunity: non-uniform selection p can lower the sampling-weighted dimension by prioritizing FO-dominant updates, at the cost of a representation penalty (ζ2N‖p−1/N‖2)O(ζ^2N\|p-1/N\|^2). This principled trade-off between directional variance reduction and data representation directly drives the algorithmic design of our dimension-aware sampling policy. Problem formulation To design a principled sampling policy, we isolate the sampling-dependent terms in our multi-step trajectory descent bound (Proposition C.8). Suppressing the round index, let Ui≔(w−wiE)/hU_i (w-w_i^E)/h denote client i’s normalized local update with effective step size h≔ηEh η E, ai≔[Ui]a_i [U_i] be its expected trajectory, vi≔‖Ui−ai‖2v_i \|U_i-a_i\|^2 be its local variance, a¯≔1N∑iai a 1N _ia_i be the unweighted cohort mean, and g≔∇f(w)g ∇ f(w) represent the true global gradient. Over the constrained simplex ≔p∈ℝN:∑ipi=1,ϵ/N≤pi≤1/KC \p ^N: _ip_i=1,\;ε/N≤ p_i≤ 1/K\, minimizing the descent bound yields the ideal statistical quadratic program (QP): minp∈‖∑i=1Npiai−g‖2+LhK∑i=1Npi[vi+2‖ai−a¯‖2]. _p \; \| _i=1^Np_ia_i-g \|^2+ LhK _i=1^Np_i [v_i+2\|a_i- a\|^2 ]. (13) Problem (13) serves as an optimal reference: the first term minimizes representation bias relative to g, while the second penalizes cohort variance induced by local noise viv_i and client drift dispersion ‖ai−a¯‖2\|a_i- a\|^2. However, solving (13) online is intractable because trajectory moments ai,vi\a_i,v_i\ and gradient g are unobservable prior to local execution. Scalar relaxation and water-filling structure. To decouple directional cross-terms and eliminate the unknown gradient g, we relax the squared bias via Cauchy–Schwarz: ‖∑ipiai−g‖2≤2NHa‖p−u‖2+2‖a¯−g‖2\| _ip_ia_i-g\|^2≤ 2NH_a\|p-u\|^2+2\| a-g\|^2, where Ha≔1N∑i‖ai−a¯‖2H_a 1N _i\|a_i- a\|^2 measures empirical update dispersion and u≔[1/N,…,1/N]⊤u [1/N,…,1/N] . Dropping the p-independent offset 2‖a¯−g‖22\| a-g\|^2 transforms (13) into a separable scalar program: minp∈A‖p−u‖2+c⊤p,whereA≔2NHa,ci≔LhK[vi+2‖ai−a¯‖2]. _p \;A\|p-u\|^2+c p, A 2NH_a, c_i LhK [v_i+2\|a_i- a\|^2 ]. (14) When A>0A>0, the KKT conditions admit a closed-form water-filling solution: pi⋆=clip[ϵ/N, 1/K](1N−ci+λ2A),p_i =clip_[ε/N,\,1/K] ( 1N- c_i+λ2A ), (15) where scalar dual multiplier λ enforces ∑ipi⋆=1 _ip_i =1. Appendix D.1 provides the full relaxation and KKT derivations. Equation (15) reveals the qualitative behavior of optimal sampling: it prioritizes clients with smaller noise costs cic_i up to the capacity cap 1/K1/K, while pulling coordinates toward uniform 1/N1/N to bound representation error. The practical HO-FL policy. Evaluating runtime costs cic_i remains costly as it requires dynamic variance tracking. To achieve zero runtime overhead while retaining the water-filling principle, we substitute cic_i with the static dimension score si≔diZO/qs_i d_i^ZO/q, grounded in the (diZO+1)/q(d_i^ZO+1)/q variance penalty in Lemma C.5. Under this proxy, minimizing noise alone (minp∈s⊤p _p s p) is greedily solved by saturating the K clients with the smallest scores at 1/K1/K, yielding ri≔1Ki∈Rr_i 1K1\i∈ R\, where R denotes the greedy subset. Directly parameterizing the trade-off line segment connecting the representation-optimal distribution u and the noise-optimal distribution r produces our closed-form HO-FL policy: pi=(1−β)1N+βri,β∈[0,1).p_i=(1-β) 1N+β r_i, β∈[0,1). (16) Equation (16) automatically satisfies the bounded simplex constraints with pi≤1/Kp_i≤ 1/K and coverage floor pi≥(1−β)/Np_i≥(1-β)/N. Here, hyperparameter β serves as an explicit control knob governing the theoretical trade-off in Theorem 4.1: setting β=0β=0 recovers uniform sampling and eliminates representation bias (‖p−u‖2=0\|p-u\|^2=0), whereas increasing β tilts selection toward FO-dominant clients, actively suppressing directional variance while analytically bounding representation bias by N‖p−u‖2=β2(N/K−1)N\|p-u\|^2=β^2(N/K-1). Practical guidance for selecting a fixed β provided in Appendix D.2. 6 Experiments In this section, we empirically evaluate HO-FL across diverse foundation models and downstream tasks, centered on three core research questions: • RQ1 (Task Performance): Can HO-FL close the accuracy gap to full first-order federated optimization while accommodating heterogeneous client order boundary? • RQ2 (Memory Feasibility): How effectively does hybrid-order local training reduce client peak activation memory, and what are the system-wide resource footprints? • RQ3 (Heterogeneity & Trade-offs): How does the dimension-aware sampling policy (β) interact with data non-IIDness (α), and does it validate our theoretical convergence bound? 6.1 Experimental Setup Federated environment and baselines. We simulate an edge system of N=30N=30 clients with cohort size K=6K=6 over 160 communication rounds. Each active client executes E=5E=5 local AdamW steps. To emulate memory heterogeneity, clients are evenly partitioned across three label-independent boundary tiers (Table 1). We benchmark HO-FL (β=0.75β=0.75) against full-FO (FedAvg (McMahan et al., 2017), FedProx (Li et al., 2020)), pure-ZO DeComFL (Li et al., 2025), and the uniform sampling version of HO-FL (β=0β=0). Implementation details are in Appendix E. Table 1: Boundary tiers. Order Boundary setting across 3 tiers. Model Layers Order Boundary T1 T2 T3 OPT-125M 12 3 6 9 Qwen2.5-1.5B 28 8 16 24 SmolLM3-3B 36 9 18 27 Models and benchmarks. We evaluate HO-FL by fine-tuning OPT-125M, Qwen2.5-1.5B, and SmolLM3-3B (Zhang et al., 2022; Qwen et al., 2024; Hugging Face, 2025) with LoRA adapters (Hu et al., 2022) . The benchmarks include SST-2, BoolQ, SciQ, SQuAD 1.1, and AG News (Socher et al., 2013; Clark et al., 2019; Welbl et al., 2017; Rajpurkar et al., 2016; Zhang et al., 2015). Tasks are evaluated under IID client partitions, while AG News additionally incorporates Dirichlet label-skewed partitions Dir(α)Dir(α) with α∈1,0.5,0.1α∈\1,0.5,0.1\ to evaluate non-IID robustness. 6.2 Empirical Evaluation and Analysis Table 2: Language-task accuracy (%). Reported as mean(SD)mean_(SD) across three seeds. All tasks are IID except AG News (α=0.5α=0.5). Model Task FedAvg FedProx DeComFL HO-FL (β=0β=0) HO-FL (β=0.75β=0.75) OPT-125M SST-2 87.04(0.11)87.04_(0.11) 87.12(0.13)87.12_(0.13) 51.95(0.00)51.95_(0.00) 86.01(0.11)86.01_(0.11) 86.54(0.24)86.54_(0.24) BoolQ 59.84(0.37)59.84_(0.37) 59.85(0.46)59.85_(0.46) 55.87(0.05)55.87_(0.05) 57.91(0.18)57.91_(0.18) 59.44(0.96)59.44_(0.96) SciQ 25.50(1.30)25.50_(1.30) 25.43(1.04)25.43_(1.04) 25.53(0.85)25.53_(0.85) 24.97(0.80)24.97_(0.80) 26.10(1.64)26.10_(1.64) AG News (α=0.5α=0.5) 78.24(2.97)78.24_(2.97) 78.29(2.99)78.29_(2.99) 23.38(0.08)23.38_(0.08) 59.58(2.24)59.58_(2.24) 62.85(6.52)62.85_(6.52) Qwen2.5-1.5B SST-2 93.27(0.26)93.27_(0.26) 93.43(0.26)93.43_(0.26) 71.25(0.18)71.25_(0.18) 91.09(0.18)91.09_(0.18) 92.81(0.07)92.81_(0.07) BoolQ 79.06(0.19)79.06_(0.19) 79.11(0.08)79.11_(0.08) 72.44(0.12)72.44_(0.12) 77.00(0.40)77.00_(0.40) 78.72(0.24)78.72_(0.24) SciQ 90.83(0.32)90.83_(0.32) 90.87(0.31)90.87_(0.31) 88.87(0.46)88.87_(0.46) 90.17(0.25)90.17_(0.25) 90.47(0.29)90.47_(0.29) AG News (α=0.5α=0.5) 88.15(0.21)88.15_(0.21) 88.23(0.17)88.23_(0.17) 55.83(0.08)55.83_(0.08) 85.51(0.43)85.51_(0.43) 86.69(0.56)86.69_(0.56) SmolLM3-3B SST-2 94.04(0.11)94.04_(0.11) 94.07(0.07)94.07_(0.07) 63.38(0.18)63.38_(0.18) 93.16(0.07)93.16_(0.07) 94.04(0.30)94.04_(0.30) BoolQ 86.15(0.05)86.15_(0.05) 86.01(0.14)86.01_(0.14) 82.59(0.07)82.59_(0.07) 84.32(0.12)84.32_(0.12) 85.37(0.43)85.37_(0.43) SciQ 90.50(0.17)90.50_(0.17) 90.57(0.35)90.57_(0.35) 89.13(0.38)89.13_(0.38) 89.57(0.35)89.57_(0.35) 89.77(0.21)89.77_(0.21) AG News (α=0.5α=0.5) 89.38(0.14)89.38_(0.14) 89.37(0.18)89.37_(0.18) 74.70(0.05)74.70_(0.05) 87.11(0.22)87.11_(0.22) 87.52(0.55)87.52_(0.55) Downstream task performance (RQ1). As shown in Table 2, HO-FL with β=0.75β=0.75 closely approaches FO baselines accuracy across language understanding tasks, remaining within 0.5%–0.8% of FedAvg on most benchmarks while surpassing DeComFL by up to 34.59% (86.54% vs. 51.95% on OPT-125M/SST-2). Moreover, dimension-aware sampling consistently outperforms uniform hybrid selection with β=0β=0, yielding up to a 3.27% accuracy gain under label skew on AG News at α=0.5α=0.5. On generative QA evaluated on SQuAD 1.1 (Table 3), where pure ZO degrades severely as DeComFL achieves only 8.54% F1 on OPT-125M, HO-FL reliably restores convergence by lifting to 52.43%. Specifically, dimension-aware sampling improves F1 over uniform hybrid selection by 1.3%–5.2% across architectures, narrowing the performance gap with FedAvg to within 0.7%–2.3%. Table 3: SQuAD 1.1 performance (%). Full-validation EM/F1, reported as mean(SD)mean_(SD) across three seeds. Model Metric FedAvg FedProx DeComFL HO-FL (β=0β=0) HO-FL (β=0.75β=0.75) OPT-125M EM 43.95(0.04)43.95_(0.04) 43.92(0.10)43.92_(0.10) 1.29(0.05)1.29_(0.05) 37.16(0.26)37.16_(0.26) 41.87(0.32)41.87_(0.32) F1 54.68(0.02)54.68_(0.02) 54.65(0.08)54.65_(0.08) 8.54(0.09)8.54_(0.09) 47.24(0.21)47.24_(0.21) 52.43(0.20)52.43_(0.20) Qwen2.5-1.5B EM 79.65(0.14)79.65_(0.14) 79.70(0.08)79.70_(0.08) 61.70(0.11)61.70_(0.11) 76.67(0.05)76.67_(0.05) 78.71(0.19)78.71_(0.19) F1 87.59(0.07)87.59_(0.07) 87.61(0.04)87.61_(0.04) 73.14(0.16)73.14_(0.16) 85.34(0.03)85.34_(0.03) 86.92(0.10)86.92_(0.10) SmolLM3-3B EM 81.52(0.10)81.52_(0.10) 81.47(0.20)81.47_(0.20) 47.98(0.09)47.98_(0.09) 77.80(0.11)77.80_(0.11) 79.54(0.20)79.54_(0.20) F1 89.93(0.10)89.93_(0.10) 89.90(0.16)89.90_(0.16) 68.43(0.04)68.43_(0.04) 87.25(0.04)87.25_(0.04) 88.57(0.13)88.57_(0.13) Figure 2: Single-client peak and system average memory footprints on Qwen2.5-1.5B. Left: measured single-client peak allocated memory across different order boundaries during local training (batch size 16 and length 1024). Right: round-averaged cohort memory demand C¯(L) C(L) across participating clients (N=30,K=6N=30,K=6) under varying β. On-device memory footprint and feasibility (RQ2). As illustrated in Figure 2, HO-FL systematically resolves the activation memory bottleneck during local training. At sequence length 1024, full BP (FedAvg) incurs a prohibitive peak memory of 40.26 GiB on Qwen2.5-1.5B. By setting the boundary at block 24 to truncate the backward computation graph, HO-FL slashes peak allocated memory to 9.08 GiB—achieving a 77.4% reduction that brings 1.5B LLM adaptation within reach of commodity 12 GiB devices. To evaluate distributed system overhead, we profile the expected cohort memory demand C¯(L)=T−1∑t∑i∈Stmi(L) C(L)=T^-1 _t _i∈ S_tm_i(L) at sequence length 1024 (Figure 2). Uniform hybrid sampling (β=0β=0) and dimension-aware sampling (β=0.75β=0.75) require an average of 117.32 GiB and 163.78 GiB per round, respectively, yielding substantial savings over homogeneous FedAvg (241.57 GiB). Although a higher β moderately increases aggregate resource utilization by prioritizing FO-dominant clients, it accelerates global convergence while strictly respecting each client’s physical memory budget. 6.3 Sampling Preference and Data Heterogeneity (RQ3) Table 4: Sampling preference and data heterogeneity. Qwen/AG News final test accuracy (%), reported as mean(SD)mean_(SD) across three seeds. Distribution β=0β=0 β=0.25β=0.25 β=0.5β=0.5 β=0.75β=0.75 IID 86.75(0.13)86.75_(0.13) 87.33(0.03)87.33_(0.03) 87.58(0.08)87.58_(0.08) 87.91(0.03)87.91_(0.03) α=1α=1 86.22(0.47)86.22_(0.47) 86.64(0.23)86.64_(0.23) 87.11(0.15)87.11_(0.15) 87.26(0.29)87.26_(0.29) α=0.5α=0.5 85.51(0.43)85.51_(0.43) 86.31(0.48)86.31_(0.48) 86.58(0.44)86.58_(0.44) 86.69(0.56)86.69_(0.56) α=0.1α=0.1 83.39(1.46)83.39_(1.46) 83.83(0.92)83.83_(0.92) 82.53(1.81)82.53_(1.81) 82.66(3.07)82.66_(3.07) On Qwen/AG News (Table 4), β=0.75β=0.75 improves over HO-FL (β=0β=0) by 1.17 points under IID, 1.04 points at α=1α=1, and 1.18 points at α=0.5α=0.5, with gains in every seed. Under stronger skew (α=0.1α=0.1), 0.250.25 is preferable: 0.50.5 and 0.750.75 fall below HO-FL (β=0β=0) by 0.87 and 0.74 points, and 0.750.75 has SD 3.07. This supports the trade-off between update quality and data representation as shown in Theorem 4.1. Appendix E.6 further reports E/q/μE/q/μ ablations on OPT/SST-2. 7 Conclusion In this paper, we presented HO-FL, the first fully federated hybrid-order learning framework designed for heterogeneous edge devices. By executing ZO optimization on the bottom segment and FO backpropagation on the top segment entirely on each client, HO-FL enables memory-constrained devices to adapt their local order boundaries to memory budgets. Our theoretical analysis characterizes the convergence behavior of HO-FL and introduces a principled dimension-aware client sampling policy that balances directional gradient variance against global data representation. Extensive empirical evaluations across diverse tasks demonstrate that HO-FL achieves performance competitive with full first-order federated learning while slashing client peak memory by up to 77.4%. Promising avenues for future work include extending HO-FL to support runtime-adaptive order boundaries under fluctuating edge resources and developing dynamic schedules for β . References Chen et al. (2026) Q. Chen, X. Wu, Y. Wang, and X. Chen HO-SFL: hybrid-order split federated learning with backprop-free clients and dimension-free aggregation. In International Conference on Machine Learning, External Links: Link Cited by: 2nd item, Table 5, §C.1, §1, §1, §2, §3, §4.1, §4.2. Clark et al. (2019) C. Clark, K. Lee, M. Chang, T. Kwiatkowski, M. Collins, and K. Toutanova BoolQ: exploring the surprising difficulty of natural yes/no questions. In Proceedings of the Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, p. 2924–2936. External Links: Document, Link Cited by: 2nd item, §6.1. Fang et al. (2022) W. Fang, Z. Yu, Y. Jiang, Y. Shi, C. N. Jones, and Y. Zhou Communication-efficient stochastic zeroth-order optimization for federated learning. IEEE Transactions on Signal Processing 70, p. 5058–5073. External Links: Document Cited by: §1, §2, §4.2. Gandhi et al. (2006) R. Gandhi, S. Khuller, S. Parthasarathy, and A. Srinivasan Dependent rounding and its applications to approximation algorithms. Journal of the ACM 53 (3), p. 324–360. External Links: Document Cited by: §3. Hu et al. (2022) E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen LoRA: low-rank adaptation of large language models. In International Conference on Learning Representations, External Links: Link Cited by: §E.2, §6.1. Hugging Face (2025) Hugging Face SmolLM3-3B model card. External Links: Link Cited by: §6.1. Karimireddy et al. (2020) S. P. Karimireddy, S. Kale, M. Mohri, S. Reddi, S. Stich, and A. T. Suresh SCAFFOLD: stochastic controlled averaging for federated learning. In Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 119, p. 5132–5143. External Links: Link Cited by: §C.1, §4.1. Khodak et al. (2021) M. Khodak, R. Tu, T. Li, L. Li, M. Balcan, V. Smith, and A. Talwalkar Federated hyperparameter tuning: challenges, baselines, and connections to weight-sharing. In Advances in Neural Information Processing Systems, Vol. 34. External Links: Link Cited by: §D.2. Kou et al. (2026) Z. Kou, Z. Chen, J. Yang, and C. Shen Lean clients, full accuracy: hybrid zeroth- and first-order split federated learning. arXiv preprint arXiv:2601.09076. External Links: Link Cited by: 3rd item, Table 5, §1, §2. Li et al. (2022) Q. Li, Y. Diao, Q. Chen, and B. He Federated learning on Non-IID data silos: an experimental study. In 2022 IEEE 38th International Conference on Data Engineering (ICDE), p. 965–978. External Links: Document, Link Cited by: 4th item. Li et al. (2020) T. Li, A. K. Sahu, M. Zaheer, M. Sanjabi, A. Talwalkar, and V. Smith Federated optimization in heterogeneous networks. Proceedings of Machine Learning and Systems 2, p. 429–450. External Links: Link Cited by: 2nd item, Table 1. Li et al. (2025) Z. Li, B. Ying, Z. Liu, C. Dong, and H. Yang Achieving dimension-free communication in federated learning via zeroth-order optimization. In International Conference on Learning Representations, External Links: Link Cited by: 3rd item, §1, §2, §4.2, Table 1. Li et al. (2026) Z. Li, B. Ying, Z. Liu, C. Dong, and H. Yang Converge faster, talk less: hessian-informed federated zeroth-order optimization. In International Conference on Learning Representations, External Links: Link Cited by: §1, §2. Lnu et al. (2026) A. Lnu, Z. Li, D. Liang, C. Huang, R. Li, and H. Yang HOSL: hybrid-order split learning for memory-constrained edge training. In International Symposium on Modeling and Optimization in Mobile, Ad Hoc, and Wireless Networks, p. 1–8. External Links: Document, Link Cited by: 1st item, Table 5, §1, §2. Loshchilov and Hutter (2019) I. Loshchilov and F. Hutter Decoupled weight decay regularization. In International Conference on Learning Representations, External Links: Link Cited by: §B.3, 2nd item. Malladi et al. (2023) S. Malladi, T. Gao, E. Nichani, A. Damian, J. D. Lee, D. Chen, and S. Arora Fine-tuning language models with just forward passes. In Advances in Neural Information Processing Systems, Vol. 36, p. 53038–53075. External Links: Link Cited by: §2. McMahan et al. (2017) B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. Agüera y Arcas Communication-efficient learning of deep networks from decentralized data. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 54, p. 1273–1282. External Links: Link Cited by: 1st item, §1, Table 1. Qwen et al. (2024) Qwen, A. Yang, B. Yang, B. Zhang, B. Hui, B. Zheng, B. Yu, C. Li, D. Liu, F. Huang, et al. Qwen2.5 technical report. arXiv preprint arXiv:2412.15115. External Links: Link Cited by: §6.1. Rajpurkar et al. (2016) P. Rajpurkar, J. Zhang, K. Lopyrev, and P. Liang SQuAD: 100,000+ questions for machine comprehension of text. In Proceedings of the Conference on Empirical Methods in Natural Language Processing, p. 2383–2392. External Links: Document, Link Cited by: 5th item, §6.1. Socher et al. (2013) R. Socher, A. Perelygin, J. Wu, J. Chuang, C. D. Manning, A. Y. Ng, and C. Potts Recursive deep models for semantic compositionality over a sentiment treebank. In Proceedings of the Conference on Empirical Methods in Natural Language Processing, p. 1631–1642. External Links: Link Cited by: 1st item, §6.1. Wang et al. (2024) J. Wang, R. Das, G. Joshi, S. Kale, Z. Xu, and T. Zhang On the unreasonable effectiveness of federated averaging with heterogeneous data. Transactions on Machine Learning Research. External Links: Link Cited by: §C.3, §4.2. Welbl et al. (2017) J. Welbl, N. F. Liu, and M. Gardner Crowdsourcing multiple choice science questions. In Proceedings of the 3rd Workshop on Noisy User-generated Text, L. Derczynski, W. Xu, A. Ritter, and T. Baldwin (Eds.), Copenhagen, Denmark, p. 94–106. External Links: Link, Document Cited by: 3rd item, §6.1. Wu et al. (2026) Y. Wu, C. Tian, J. Li, H. Sun, K. Tam, Z. Zhou, H. Liao, J. Xiong, Z. Guo, L. Li, and C. Xu A survey on federated fine-tuning of large language models. Transactions on Machine Learning Research. External Links: Link Cited by: §1. Zhang et al. (2022) S. Zhang, S. Roller, N. Goyal, M. Artetxe, M. Chen, S. Chen, C. Dewan, M. Diab, X. Li, X. V. Lin, et al. OPT: open pre-trained transformer language models. arXiv preprint arXiv:2205.01068. External Links: Link Cited by: §6.1. Zhang et al. (2015) X. Zhang, J. Zhao, and Y. LeCun Character-level convolutional networks for text classification. In Advances in Neural Information Processing Systems, Vol. 28, p. 649–657. External Links: Link Cited by: 4th item, §6.1. Appendix Contents Appendix A A Unified Formulation of Hybrid-Order Optimization This section formalizes the generalized hybrid-order (HO) optimization framework referenced in Section 2, unifying existing distributed hybrid-order paradigms under a common structural formulation. A.1 Structural Decomposition and Forward Propagation Consider a deep neural network partitioned across a specified cut layer into a bottom sub-network fZO(⋅,wZO)f_ZO(·;w^ZO) and a top sub-network fFO(⋅,wFO)f_FO(·;w^FO). The parameter space is canonically decoupled as: w w =[(wZO)⊤,(wFO)⊤]⊤∈ℝd,dZO+dFO=d, = [(w^ZO) ,(w^FO) ] ^d, d^ZO+d^FO=d, (17) wZO w^ZO ∈ℝdZO,wFO∈ℝdFO. ^d^ZO, w^FO ^d^FO. Given an input sample x with target label y, forward execution proceeds sequentially through the boundary activation z to produce the task prediction y y: z z =fZO(x;wZO),y^=fFO(z;wFO), =f_ZO(x;w^ZO), y=f_FO(z;w^FO), (18) ℓ(w,x,y) (w;x,y) =ℓ(y^,y). = ( y,y). In all hybrid architectures, the top segment wFOw^FO is optimized via standard first-order backpropagation, yielding the exact gradient gFO=∇wFOℓ(y^,y)g^FO= _w^FO ( y,y). The bottom segment wZOw^ZO is updated via zeroth-order perturbations evaluated over a local surrogate objective ℒlocal(v)L_local(v) at candidate parameters v. A.2 Taxonomy of Local Surrogate Objectives The fundamental distinction among existing distributed hybrid-order methods lies in how the surrogate objective ℒlocal(v)L_local(v) is constructed: • Exact End-to-End Evaluation (HOSL (Lnu et al., 2026)): HOSL perturbs the bottom segment while evaluating the entire forward graph: ℒlocalHOSL(v)=ℓ(fFO(fZO(x,v),wFO),y).L_local^HOSL(v)= (f_FO(f_ZO(x;v);w^FO),y ). (19) Each perturbation query requires transmitting the perturbed activation fZO(x,v)f_ZO(x;v) to the server and returning the scalar loss, incurring substantial per-iteration communication. • First-Order Activation Feedback (HO-SFL (Chen et al., 2026)): HO-SFL freezes the cut-layer activation gradient δ=∇zℓ(fFO(z,wFO),y)δ= _z (f_FO(z;w^FO),y) evaluated at the unperturbed anchor, forming a linear inner-product surrogate: ℒlocalHO-SFL(v)=⟨δ,fZO(x,v)⟩.L_local^HO -SFL(v)= δ,f_ZO(x;v) . (20) This formulation strictly confines all perturbation forward passes to the bottom sub-network, eliminating activation communication during the ZO search phase. • Auxiliary Network Decoupling (HERON-SFL (Kou et al., 2026)): HERON-SFL attaches a lightweight client-side auxiliary classifier h(⋅,θ)h(·;θ) to approximate the top segment’s guidance locally: ℒlocalHERON(v,θ)=ℓ(h(fZO(x,v),θ),y).L_local^HERON(v,θ)= (h(f_ZO(x;v);θ),y ). (21) Both v and the auxiliary parameters θ are updated locally via ZO perturbations to minimize cross-node dependency, at the cost of auxiliary architectural modeling bias. A.3 Architectural Comparison: Split Computing vs. HO-FL As summarized in Table 5, HOSL, HO-SFL, and HERON-SFL use Split Learning (SL) or Split Federated Learning (SFL) workflows, with cross-node activation transmissions between clients and the server. In contrast, HO-FL adopts the computational efficiency of the linear surrogate (20) but breaks away from the split paradigm: it executes the complete hybrid-order pipeline on-device, natively accommodates client-adaptive boundaries bib_i matched to heterogeneous physical memory constraints, and aggregates model updates after E local training steps using seed replay for the ZO segment and transmitted parameter differences for the FO segment. Table 5: Structural comparison of distributed hybrid-order paradigms. Method System Paradigm Boundary Flexibility Activation Uploads Local Steps per Aggregation HOSL (Lnu et al., 2026) Split Learning Fixed (Single Client) 2q+12q+1 uploads / step N/A (Single Client) HO-SFL (Chen et al., 2026) Split FedLearning Identical 1 upload / step 1 HERON-SFL (Kou et al., 2026) Split FedLearning Identical Every k local steps h≥1h≥ 1 HO-FL (Ours) Federated Learning Heterogeneous (Adaptive) None (Fully On-Device) ≥ E≥ 1 Appendix B Algorithm and Communication Details B.1 Fixed-size dependent rounding Given normalized inclusion probabilities pi∈[0,1/K]p_i∈[0,1/K] with ∑ipi=1 _ip_i=1, Algorithm 2 rounds x=Kpx=Kp to a binary vector. It preserves [Xi]=xiE[X_i]=x_i and ∑iXi=K _iX_i=K exactly and satisfies [XiXj]≤xixjE[X_iX_j]≤ x_ix_j for i≠ji≠ j. One randomized permutation is drawn before the pairwise procedure; the next two fractional entries are taken in this order. When all inputs are equal, permutation symmetry makes the output uniform over all K-subsets. Algorithm 2 Function: DepRound 1: x∈[0,1]Nx∈[0,1]^N, ∑ixi=K∈ℕ _ix_i=K . 2: function DepRound(x) 3: π←π← a uniformly random permutation of 1,…,N\1,…,N\. 4: while |i:0<xi<1|≥2|\i:0<x_i<1\|≥ 2 do 5: i,j←i,j← the first two fractional indices in order π. 6: δ+←min(1−xi,xj) _+← (1-x_i,x_j); δ−←min(xi,1−xj) _-← (x_i,1-x_j). 7: a∼Uniform(0,1)a (0,1). 8: if a<δ−/(δ++δ−)a< _-/( _++ _-) then 9: (xi,xj)←(xi+δ+,xj−δ+)(x_i,x_j)←(x_i+ _+,x_j- _+). 10: else 11: (xi,xj)←(xi−δ−,xj+δ−)(x_i,x_j)←(x_i- _-,x_j+ _-). 12: end if 13: end while 14: return S=i:xi=1S=\i:x_i=1\. ⊳ Used by Algorithm 1 15: end function Every step preserves the sum, preserves each coordinate in expectation, and makes at least one coordinate integral. For a pair updated together, its product has conditional expected decrement δ+δ− _+ _-; when only one coordinate of a fixed pair is updated, that pair’s product is a martingale. Thus pair products are supermartingales, giving the stated nonpositive inclusion covariances. Because the sum is an integer, one fractional coordinate cannot remain. Implementations snap values within floating-point tolerance of 00 or 11 and check the final cardinality. The probability rule used in the main experiments. Let R contain the K clients with the smallest diZO/qd_i^ZO/q, with ties resolved by a seeded random permutation drawn once at run start and retained for stable sorting throughout the run. For ui=1/Nu_i=1/N, ri=i∈RK,pi=(1−β)ui+βri,St=DepRound(Kp).r_i= 1\i∈ R\K, p_i=(1-β)u_i+β r_i, S_t=DepRound(Kp). (22) If all resource scores are equal, set r=ur=u. The IID main comparisons use β=0.75β=0.75; the sampling sensitivity study uses β∈0,0.25,0.5,0.75β∈\0,0.25,0.5,0.75\, and the OPT/SST-2 local-training ablations fix β=0.5β=0.5. The resource ranking and p are fixed before training and held constant across rounds. Dependent rounding uses a fresh random permutation each round to draw the participants. The resulting p always satisfies pi≤1/Kp_i≤ 1/K and pi≥(1−β)/Np_i≥(1-β)/N. B.2 Local updates and server reconstruction Algorithm 3 implements the E-step local update in (7) and returns the directional scalars and top-segment parameter difference. Transient quantities omit round and step indices and are recomputed at each step; ℓ(y^i,y) ( y_i,y) denotes the prediction-space form of ℓ(wi,ξ) (w_i;ξ). Direction streams are independent across clients and steps and separate from minibatch and random-layer streams. For full-FO clients, the ZO loop and scalar record are empty. Algorithm 3 Function: ClientUpdate 1: bi,E,q,μ,ηb_i,E,q,μ,η and dispatched seedi,t,jee,j\seed_i,t,j^e\_e,j from Algorithm 1. 2: function ClientUpdate(i,t,wti,t,w_t) 3: wi,t0←wtw_i,t^0← w_t. 4: for e=0,…,E−1e=0,…,E-1 do 5: ξi,te=(x,y)∼i _i,t^e=(x,y) _i. 6: // Forward and top-segment BP 7: zi←fi,ZO(x,wi,tZO,e)z_i← f_i,ZO(x;w_i,t^ZO,e). ⊳ Forward without caching activation 8: y^i←fi,FO(zi,wi,tFO,e) y_i← f_i,FO(z_i;w_i,t^FO,e). 9: giFO←∇wi,tFO,eℓ(y^i,y)g_i^FO← _w_i,t^FO,e ( y_i,y), δi←∇ziℓ(y^i,y) _i← _z_i ( y_i,y). ⊳ Eq. (3) 10: // Surrogate and directional scalars 11: ℒi(v)←⟨δi,fi,ZO(x,v)⟩L_i(v)← _i,f_i,ZO(x;v) . ⊳ Eq. (4) 12: g^iZO←0 g_i^ZO← 0. 13: if diZO>0d_i^ZO>0 then 14: for j=1,…,qj=1,…,q do 15: ui,j←RNG(seedi,t,je)u_i,j (seed_i,t,j^e). ⊳ Eq. (5) 16: si,t,je←ℒi(wi,tZO,e+μui,j)−ℒi(wi,tZO,e)μs_i,t,j^e← L_i(w_i,t^ZO,e+μ u_i,j)-L_i(w_i,t^ZO,e)μ. ⊳ Eq. (6) 17: g^iZO←g^iZO+1qsi,t,jeui,j g_i^ZO← g_i^ZO+ 1qs_i,t,j^eu_i,j. 18: end for 19: end if 20: // Update both segments at the same anchor 21: wi,tZO,e+1←wi,tZO,e−ηg^iZOw_i,t^ZO,e+1← w_i,t^ZO,e-η g_i^ZO. ⊳ Eq. (7) 22: wi,tFO,e+1←wi,tFO,e−ηgiFOw_i,t^FO,e+1← w_i,t^FO,e-η g_i^FO. 23: end for 24: Δi,tFO←wi,tFO,E−wi,tFO,0 _i,t^FO← w_i,t^FO,E-w_i,t^FO,0. ⊳ Eq. (8) 25: return (si,t,jee,j,Δi,tFO)(\s_i,t,j^e\_e,j, _i,t^FO). 26: end function Algorithm 4 regenerates the directions from the shared seeds and recovers the bottom update using (9), allowing Algorithm 1 to average updates across heterogeneous boundaries via (10). Algorithm 4 Function: Reconstruct 1: bi,E,q,ηb_i,E,q,η and known seedi,t,jee,j\seed_i,t,j^e\_e,j from Algorithm 1. 2: function Reconstruct(i,t,si,t,jee,j,Δi,tFOi,t,\s_i,t,j^e\_e,j, _i,t^FO) 3: ui,t,je←RNG(seedi,t,je)u_i,t,j^e (seed_i,t,j^e), e=0,…,E−1e=0,…,E-1, j=1,…,qj=1,…,q. 4: Δi,tZO←−ηq∑e=0E−1∑j=1qsi,t,jeui,t,je _i,t^ZO←- ηq _e=0^E-1 _j=1^qs_i,t,j^eu_i,t,j^e. ⊳ Eq. (9) 5: Δi,t←[Δi,tZOΔi,tFO]∈ℝd _i,t← bmatrix [rgb]0.1406,0.2891,0.4688 _i,t^ZO\\ [rgb]0.8516,0.5313,0.1953 _i,t^FO bmatrix ^d. ⊳ Shared coordinates; boundary bib_i 6: return Δi,t _i,t. 7: end function B.3 Seed Replay for Local AdamW While our theoretical analysis in Section 4 is formulated with standard SGD, practical language model fine-tuning employs adaptive optimizers such as AdamW (Loshchilov and Hutter, 2019). Under SGD, bottom-segment updates depend linearly on perturbation directions, permitting direct closed-form reconstruction via Eq. (9). In contrast, AdamW applies non-linear coordinate-wise second-moment normalization and decoupled weight decay, which prevents the server from simply summing directional scalars across local iterations. To retain dimension-free uplink communication under AdamW, HO-FL executes a deterministic seed replay at the server. Crucially, in our federated setting, each selected client’s local AdamW optimizer states (the first- and second-moment buffers m and v) are reset to zero upon selection at each round. This stateless convention eliminates the need to transmit or maintain bulky historical optimizer states across rounds, allowing the server to reconstruct the exact client trajectory using only the dispatched seeds and uploaded directional scalars: 1. Momentum Initialization: At the start of round t, the server initializes the client’s bottom parameters with the current global model segment, wi,tZO,0=wtZOw_i,t^ZO,0=w_t^ZO, and resets the moment buffers to zero (m−1=0,v−1=0m_-1=0,v_-1=0). 2. Direction Regeneration & Gradient Recovery: For each local step e∈0,…,E−1e∈\0,…,E-1\, the server regenerates the exact perturbation directions ui,t,je=RNG(seedi,t,je)u_i,t,j^e=RNG(seed_i,t,j^e) using the synchronized seed sequence, reconstructing the zeroth-order gradient estimate: g^i,tZO,e=1q∑j=1qsi,t,jeui,t,je. g_i,t^ZO,e= 1q _j=1^qs_i,t,j^eu_i,t,j^e. (23) 3. Synchronous AdamW Step: The server updates the local optimizer states and advances the parameter block via coordinate-wise AdamW: me m_e =β1me−1+(1−β1)g^i,tZO,e,ve=β2ve−1+(1−β2)(g^i,tZO,e)2, = _1m_e-1+(1- _1) g_i,t^ZO,e, v_e= _2v_e-1+(1- _2)( g_i,t^ZO,e)^2, (24) m^e m_e =me1−β1e+1,v^e=ve1−β2e+1, = m_e1- _1^e+1, v_e= v_e1- _2^e+1, (25) wi,tZO,e+1 w_i,t^ZO,e+1 =(1−ηeλwd)wi,tZO,e−ηem^ev^e+ϵopt. =(1- _e _wd)w_i,t^ZO,e- _e m_e v_e+ _opt. (26) After executing all E steps, the server obtains the final bottom parameters wi,tZO,Ew_i,t^ZO,E and derives the exact bottom update Δi,tZO=wi,tZO,E−wtZO _i,t^ZO=w_i,t^ZO,E-w_t^ZO. Appendix C Convergence Analysis This appendix provides complete, self-contained proofs for the theoretical statements in Section 4 and Section 5. C.1 Assumptions Assumption C.1 (Objective Regularity). Each local function fi:ℝd→ℝf_i:R^d is continuously differentiable and L-smooth: ‖∇fi(w)−∇fi(v)‖≤L‖w−v‖,∀w,v∈ℝd,∀i∈[N].\|∇ f_i(w)-∇ f_i(v)\|≤ L\|w-v\|, ∀ w,v ^d,\;∀ i∈[N]. (27) Furthermore, the global objective function f(w)≔1N∑i=1Nfi(w)f(w) 1N _i=1^Nf_i(w) is bounded below by finff_ . Assumption C.2 (Unbiased Reference Stochastic Gradients). For any anchor parameter w∈ℝdw ^d and a fresh stochastic minibatch ξ, conditional on the preceding local history, the reference minibatch gradient Gi(w,ξ)≔∇wℓ(w,ξ)G_i(w,ξ) _w (w;ξ) satisfies: ξ[Gi(w,ξ)] _ξ [G_i(w,ξ) ] =∇fi(w), =∇ f_i(w), (28) ξ[‖Gi(w,ξ)−∇fi(w)‖2] _ξ [\|G_i(w,ξ)-∇ f_i(w)\|^2 ] ≤σi2, ≤ _i^2, (29) ξ[‖Πi(Gi(w,ξ)−∇fi(w))‖2] _ξ [\| _i(G_i(w,ξ)-∇ f_i(w))\|^2 ] ≤σi,ZO2≤σi2, ≤ _i,ZO^2≤ _i^2, (30) where Πi≔diag(IdiZO,0diFO)∈ℝd×d _i (I_d_i^ZO,0_d_i^FO) ^d× d is the orthogonal projector onto client i’s bottom segment. Assumption C.3 (Surrogate Objective Regularity). For client i, let δi(w,ξ)≔∇ziℓ(fi,FO(zi,wFO),y) _i(w,ξ) _z_i (f_i,FO(z_i;w^FO),y) denote the boundary gradient evaluated at the unperturbed model w, where zi=fi,ZO(x,wZO)z_i=f_i,ZO(x;w^ZO). The local surrogate objective ℒi(v,w,ξ)≔⟨δi(w,ξ),fi,ZO(x,v)⟩L_i(v;w,ξ) _i(w,ξ),f_i,ZO(x;v) is LisL_i^s-smooth with respect to bottom parameters v: ‖∇vℒi(v1,w,ξ)−∇vℒi(v2,w,ξ)‖≤Lis‖v1−v2‖,∀v1,v2∈ℝdiZO.\| _vL_i(v_1;w,ξ)- _vL_i(v_2;w,ξ)\|≤ L_i^s\|v_1-v_2\|, ∀ v_1,v_2 ^d_i^ZO. (31) By the chain rule, evaluating the surrogate gradient at the unperturbed anchor recovers the exact bottom reference gradient: ∇vℒi(v,w,ξ)|v=wZO=∇wZOℓ(w,ξ). . _vL_i(v;w,ξ) |_v=w^ZO= _w^ZO (w;ξ). (32) Assumption C.4 (Affine Gradient Heterogeneity). The local data heterogeneity H(w)≔1N∑i=1N‖∇fi(w)−∇f(w)‖2H(w) 1N _i=1^N\|∇ f_i(w)-∇ f(w)\|^2 satisfies: H(w)≤aH‖∇f(w)‖2+ζ2,∀w∈ℝd.H(w)≤ a_H\|∇ f(w)\|^2+ζ^2, ∀ w ^d. (33) Discussion of assumptions. The smoothness and stochastic-gradient conditions are standard in non-convex federated optimization (Karimireddy et al., 2020; Chen et al., 2026). Assumption C.4 is equivalent to SCAFFOLD’s bounded gradient dissimilarity condition (A1), with G2=ζ2G^2=ζ^2 and B2=1+aHB^2=1+a_H. Assumption C.3 expresses the same type of boundary-gradient and network-curvature regularity used in HO-SFL (Chen et al., 2026), directly at the scalar local-objective level: a uniformly bounded boundary gradient and a Lipschitz bottom-network Jacobian imply the stated smoothness. Client-specific LisL_i^s accommodate different order boundary. Throughout, p is fixed before training and held constant across rounds. The filtration ℱtF_t contains p, wtw_t, and all history before round t sampling and local training; t[⋅]=[⋅∣ℱt]E_t[·]=E[· _t]. C.2 Moments of the Hybrid-Order Estimator Let εi(w)≔[g^i(w)∣w]−∇fi(w) _i(w) [ g_i(w) w]-∇ f_i(w) denote the estimator bias. For diZO>0d_i^ZO>0, we define the dimensional penalty and variance coefficients: κi _i ≔diZO+1q,ρi2≔μ2(Lis)2diZO, d_i^ZO+1q, _i^2 μ^2(L_i^s)^2d_i^ZO, (34) ωi2 _i^2 ≔μ2(Lis)24diZO(diZO+2)(diZO+4), μ^2(L_i^s)^24d_i^ZO(d_i^ZO+2)(d_i^ZO+4), νi _i ≔2σi2+2κiσi,ZO2+2ωi2. 2 _i^2+2 _i _i,ZO^2+2 _i^2. For full first-order clients (diZO=0d_i^ZO=0), set κi=ρi=ωi=0 _i= _i= _i=0 and νi=σi2 _i= _i^2. Lemma C.5 (Hybrid Estimator Moments). Under Assumptions C.2–C.3, the hybrid-order gradient estimator g^i(w)=[(g^iZO)⊤,(giFO)⊤]⊤ g_i(w)=[( g_i^ZO) ,(g_i^FO) ] satisfies: ‖[g^i(w)∣w]−∇fi(w)‖2 \|E[ g_i(w) w]-∇ f_i(w)\|^2 ≤ρi2, ≤ _i^2, (35) [‖g^i(w)−[g^i(w)∣w]‖2|w] [\| g_i(w)-E[ g_i(w) w]\|^2\, |\,w ] ≤2κi∥Πi∇fi(w)∥2+νi, ≤ 2 _i\| _i∇ f_i(w)\|^2+ _i, (36) [‖g^i(w)‖2|w] [\| g_i(w)\|^2\, |\,w ] ≤2∥∇fi(w)∥2+2κi∥Πi∇fi(w)∥2+νi. ≤ 2\|∇ f_i(w)\|^2+2 _i\| _i∇ f_i(w)\|^2+ _i. (37) Proof. Fix model parameter w and minibatch ξ, with independent directions u1,…,uqu_1,…,u_q drawn independently of ξ. Applying Taylor’s theorem with integral remainder to the surrogate ℒi(⋅,w,ξ)L_i(·;w,ξ) along perturbation direction u∼(0,IdiZO)u (0,I_d_i^ZO): ℒi(wZO+μu)−ℒi(wZO)μu=uu⊤∇vℒi(wZO)+R(u), L_i(w^ZO+μ u)-L_i(w^ZO)μu=u _vL_i(w^ZO)+R(u), (38) where R(u)≔∫01uu⊤[∇vℒi(wZO+τμu)−∇vℒi(wZO)]τR(u) _0^1u [ _vL_i(w^ZO+τμ u)- _vL_i(w^ZO)]dτ. By Assumption C.3, the remainder is bounded by: ‖R(u)‖≤∫01‖u‖2Lisτμ‖u‖τ=μLis2‖u‖3.\|R(u)\|≤ _0^1\|u\|^2L_i^sτμ\|u\|dτ= μ L_i^s2\|u\|^3. (39) For Gaussian perturbations u∼(0,IdiZO)u (0,I_d_i^ZO), we have [uu⊤]=IdiZOE[u ]=I_d_i^ZO and for any fixed vector x∈ℝdiZOx ^d_i^ZO: [‖1q∑j=1qujuj⊤x−x‖2]=diZO+1q‖x‖2=κi‖x‖2.E [ \| 1q _j=1^qu_ju_j x-x \|^2 ]= d_i^ZO+1q\|x\|^2= _i\|x\|^2. (40) Furthermore, by Jensen’s inequality and the standard Gaussian moment [‖u‖6]=diZO(diZO+2)(diZO+4)E[\|u\|^6]=d_i^ZO(d_i^ZO+2)(d_i^ZO+4): [‖1q∑j=1qR(uj)‖2]≤1q∑j=1q‖R(uj)‖2≤μ2(Lis)24[‖u‖6]=ωi2.E [ \| 1q _j=1^qR(u_j) \|^2 ]≤ 1q _j=1^qE\|R(u_j)\|^2≤ μ^2(L_i^s)^24E[\|u\|^6]= _i^2. (41) Decompose the estimator as g^i(w)=g^i0(w)+[(1q∑jR(uj))⊤,0⊤]⊤ g_i(w)= g_i^0(w)+[( 1q _jR(u_j)) ,0 ] , where g^i0(w) g_i^0(w) denotes the unperturbed linear component. Taking conditional expectations over the perturbation directions: u[g^i0(w)]=[(∇vℒi(wZO))⊤,(giFO)⊤]⊤=Gi(w,ξ).E_u[ g_i^0(w)]= [( _vL_i(w^ZO)) ,(g_i^FO) ] =G_i(w,ξ). (42) Taking expectations over minibatch ξ yields ξ,u[g^i0(w)]=∇fi(w)E_ξ,u[ g_i^0(w)]=∇ f_i(w). For the bias, Gaussian integration by parts gives u[ℒi(v+μu)−ℒi(v)μu]=u[∇vℒi(v+μu)]E_u [ L_i(v+μ u)-L_i(v)μu ]=E_u[ _vL_i(v+μ u)]. Therefore: ‖[g^i(w)∣w]−∇fi(w)‖ \|E[ g_i(w) w]-∇ f_i(w)\| =‖ξ,u[∇vℒi(wZO+μu)−∇vℒi(wZO)]‖ =\|E_ξ,u[ _vL_i(w^ZO+μ u)- _vL_i(w^ZO)]\| (43) ≤μLis‖u‖≤μLisdiZO=ρi, ≤μ L_i^sE\|u\|≤μ L_i^s d_i^ZO= _i, proving (35). For variance, applying ‖x+y‖2≤2‖x‖2+2‖y‖2\|x+y\|^2≤ 2\|x\|^2+2\|y\|^2 to the centered estimator: [‖g^i(w)−g^i(w)‖2] [\| g_i(w)-E g_i(w)\|^2 ] ≤2[‖g^i0(w)−∇fi(w)‖2]+2[‖1q∑j=1qR(uj)‖2] ≤ 2E [\| g_i^0(w)-∇ f_i(w)\|^2 ]+2E [ \| 1q _j=1^qR(u_j) \|^2 ] ≤2(κi‖ΠiGi(w,ξ)‖2+‖Gi(w,ξ)−∇fi(w)‖2)+2ωi2 ≤ 2 ( _iE\| _iG_i(w,ξ)\|^2+E\|G_i(w,ξ)-∇ f_i(w)\|^2 )+2 _i^2 ≤2(κi∥Πi∇fi(w)∥2+κiσi,ZO2+σi2)+2ωi2 ≤ 2 ( _i\| _i∇ f_i(w)\|^2+ _i _i,ZO^2+ _i^2 )+2 _i^2 =2κi∥Πi∇fi(w)∥2+νi, =2 _i\| _i∇ f_i(w)\|^2+ _i, (44) which establishes (36). Eq. (37) follows directly by applying ‖x+y‖2≤2‖x‖2+2‖y‖2\|x+y\|^2≤ 2\|x\|^2+2\|y\|^2 to the uncentered decomposition: ‖g^i(w)‖2 \| g_i(w)\|^2 ≤2‖g^i0(w)‖2+2‖1q∑j=1qR(uj)‖2 ≤ 2E\| g_i^0(w)\|^2+2E \| 1q _j=1^qR(u_j) \|^2 ≤2∥∇fi(w)∥2+2κi∥Πi∇fi(w)∥2+2σi2+2κiσi,ZO2+2ωi2 ≤ 2\|∇ f_i(w)\|^2+2 _i\| _i∇ f_i(w)\|^2+2 _i^2+2 _i _i,ZO^2+2 _i^2 =2∥∇fi(w)∥2+2κi∥Πi∇fi(w)∥2+νi. =2\|∇ f_i(w)\|^2+2 _i\| _i∇ f_i(w)\|^2+ _i. For full first-order clients, the reference-gradient bounds give the stated moments directly. ∎ C.3 Multi-Step Local Trajectory Analysis We analyze the local updates in (7), initialized at wi,t0=wtw_i,t^0=w_t, with a common constant step size η>0η>0 and E local steps. Each step uses a fresh minibatch and independent Gaussian directions, with independent client training streams. Let h≔ηEh η E and define the normalized trajectory and its conditional moments: Ui,t U_i,t ≔wt−wi,tEh,ai,t≔t[Ui,t], w_t-w_i,t^Eh, a_i,t _t[U_i,t], (45) vi,t v_i,t ≔t‖Ui,t−ai,t‖2. _t\|U_i,t-a_i,t\|^2. These statistics describe the local trajectory client i would execute from wtw_t, whether or not it is selected. Write a¯t≔1N∑i=1Nai,t,Bt(p)≔∑i=1Npiai,t−∇f(wt). a_t 1N _i=1^Na_i,t, B_t(p) _i=1^Np_ia_i,t-∇ f(w_t). (46) Direct model averaging gives wt+1=wt−hMtw_t+1=w_t-hM_t, where Mt=1K∑i∈StUi,tM_t= 1K _i∈ S_tU_i,t. At round t, write gi,t=∇fi(wt)g_i,t=∇ f_i(w_t) and gt=∇f(wt)g_t=∇ f(w_t). Define the local drift and average estimator bias along client i’s trajectory as di,t d_i,t ≔1E∑e=0E−1t[∇fi(wi,te)−gi,t], 1E _e=0^E-1E_t[∇ f_i(w_i,t^e)-g_i,t], (47) ε¯i,t _i,t ≔1E∑e=0E−1tεi(wi,te). 1E _e=0^E-1E_t _i(w_i,t^e). Then ai,t=gi,t+di,t+ε¯i,ta_i,t=g_i,t+d_i,t+ _i,t, and the aggregate bias separates three sources of error: Bt(p)=∑i(pi−ui)gi,t⏟selection skew+∑ipidi,t⏟local trajectory drift+∑ipiε¯i,t⏟hybrid approximation bias,B_t(p)= _i(p_i-u_i)g_i,t_selection skew+ _ip_id_i,t_local trajectory drift+ _ip_i _i,t_hybrid approximation bias, (48) where ui=1Nu_i= 1N. Retaining these vector sums preserves cancellation across clients (Wang et al., 2024). Fix round index t and suppress notation t. Let τ≔η2(E−1)2τ η^2(E-1)^2, Qi≔‖gi‖2+κi‖Πigi‖2Q_i \|g_i\|^2+ _i\| _ig_i\|^2, and define the drift radius: Di≔2τ(4Qi+νi).D_i 2τ(4Q_i+ _i). (49) Lemma C.6 (Local Trajectory Drift and Statistics). Under Assumptions C.1–C.3, if the step size satisfies 4(1+κi)L2τ≤124(1+ _i)L^2τ≤ 12, then for all e∈0,…,E−1e∈\0,…,E-1\: t‖wie−wt‖2≤Di.E_t\|w_i^e-w_t\|^2≤ D_i. (50) Furthermore, the multi-step drift vector did_i and expected update aia_i satisfy: ‖di‖2≤L2Di,‖ε¯i‖2≤ρi2,‖ai−gi‖2≤2L2Di+2ρi2,\|d_i\|^2≤ L^2D_i, \| _i\|^2≤ _i^2, \|a_i-g_i\|^2≤ 2L^2D_i+2 _i^2, (51) and the multi-step local variance satisfies: vi≤3L2Di+3ρi2+3E(4κi‖Πigi‖2+4κiL2Di+νi).v_i≤ 3L^2D_i+3 _i^2+ 3E (4 _i\| _ig_i\|^2+4 _iL^2D_i+ _i ). (52) Proof. Let Ri≔max0≤e<Et‖wie−wt‖2R_i _0≤ e<EE_t\|w_i^e-w_t\|^2. By Lemma C.5, L-smoothness, and ‖Πix‖≤‖x‖\| _ix\|≤\|x\|: t‖g^i(wie)‖2 _t\| g_i(w_i^e)\|^2 ≤2t∥∇fi(wie)∥2+2κit∥Πi∇fi(wie)∥2+νi ≤ 2E_t\|∇ f_i(w_i^e)\|^2+2 _iE_t\| _i∇ f_i(w_i^e)\|^2+ _i ≤4∥∇fi(wt)∥2+4κi∥Πi∇fi(wt)∥2+νi+4(1+κi)L2t∥wie−wt∥2 ≤ 4\|∇ f_i(w_t)\|^2+4 _i\| _i∇ f_i(w_t)\|^2+ _i+4(1+ _i)L^2E_t\|w_i^e-w_t\|^2 ≤4Qi+νi+4(1+κi)L2Ri. ≤ 4Q_i+ _i+4(1+ _i)L^2R_i. (53) Expanding the local iteration wie−wt=−η∑s=0e−1g^i(wis)w_i^e-w_t=-η _s=0^e-1 g_i(w_i^s) and applying Cauchy–Schwarz: t‖wie−wt‖2≤η2e∑s=0e−1t‖g^i(wis)‖2≤η2(E−1)2[4Qi+νi+4(1+κi)L2Ri].E_t\|w_i^e-w_t\|^2≤η^2e _s=0^e-1E_t\| g_i(w_i^s)\|^2≤η^2(E-1)^2 [4Q_i+ _i+4(1+ _i)L^2R_i ]. (54) Taking the maximum over e<Ee<E on both sides: Ri≤τ(4Qi+νi)+4(1+κi)L2τRi≤τ(4Qi+νi)+12Ri.R_i≤τ(4Q_i+ _i)+4(1+ _i)L^2τ R_i≤τ(4Q_i+ _i)+ 12R_i. (55) Subtracting 12Ri 12R_i and multiplying by 2 proves (50). Next, by L-smoothness and Jensen’s inequality: ‖di‖2=‖1E∑e=0E−1t[∇fi(wie)−gi]‖2≤1E∑e=0E−1L2t‖wie−wt‖2≤L2Di.\|d_i\|^2= \| 1E _e=0^E-1E_t[∇ f_i(w_i^e)-g_i] \|^2≤ 1E _e=0^E-1L^2E_t\|w_i^e-w_t\|^2≤ L^2D_i. (56) Since ‖εi(w)‖2≤ρi2\| _i(w)\|^2≤ _i^2 uniformly, ‖ε¯i‖2≤ρi2\| _i\|^2≤ _i^2. Using ai−gi=di+ε¯ia_i-g_i=d_i+ _i and ‖x+y‖2≤2‖x‖2+2‖y‖2\|x+y\|^2≤ 2\|x\|^2+2\|y\|^2 proves (51). To bound the variance vi=t‖Ui−ai‖2v_i=E_t\|U_i-a_i\|^2, decompose Ui=1E∑e=0E−1g^i(wie)U_i= 1E _e=0^E-1 g_i(w_i^e) into: Ui=gi+ℛi+i+ℳi,U_i=g_i+R_i+Z_i+M_i, (57) where ℛi≔1E∑e=0E−1(∇fi(wie)−gi)R_i 1E _e=0^E-1(∇ f_i(w_i^e)-g_i), i≔1E∑e=0E−1εi(wie)Z_i 1E _e=0^E-1 _i(w_i^e), and ℳi≔1E∑e=0E−1ℰieM_i 1E _e=0^E-1E_i^e with martingale difference (where ℱieF_i^e contains ℱtF_t and the client history before drawing step e randomness): ℰie≔g^i(wie)−[g^i(wie)∣ℱie].E_i^e g_i(w_i^e)-E[ g_i(w_i^e) _i^e]. (58) Because ℰieE_i^e is a martingale difference sequence (t[⟨ℰie,ℰis⟩]=0E_t[ _i^e,E_i^s ]=0 for e≠se≠ s): t‖ℳi‖2 _t\|M_i\|^2 =1E2∑e=0E−1t‖ℰie‖2 = 1E^2 _e=0^E-1E_t\|E_i^e\|^2 ≤1E2∑e=0E−1(2κit∥Πi∇fi(wie)∥2+νi) ≤ 1E^2 _e=0^E-1 (2 _iE_t\| _i∇ f_i(w_i^e)\|^2+ _i ) ≤1E(4κi‖Πigi‖2+4κiL2Di+νi). ≤ 1E (4 _i\| _ig_i\|^2+4 _iL^2D_i+ _i ). (59) Since ai=gi+tℛi+tia_i=g_i+E_tR_i+E_tZ_i, centering yields: Ui−ai=(ℛi−tℛi)+(i−ti)+ℳiU_i-a_i=(R_i-E_tR_i)+(Z_i-E_tZ_i)+M_i. Applying ‖x+y+z‖2≤3(‖x‖2+‖y‖2+‖z‖2)\|x+y+z\|^2≤ 3(\|x\|^2+\|y\|^2+\|z\|^2): vi≤3t‖ℛi‖2+3t‖i‖2+3t‖ℳi‖2≤3L2Di+3ρi2+3E(4κi‖Πigi‖2+4κiL2Di+νi),v_i≤ 3E_t\|R_i\|^2+3E_t\|Z_i\|^2+3E_t\|M_i\|^2≤ 3L^2D_i+3 _i^2+ 3E (4 _i\| _ig_i\|^2+4 _iL^2D_i+ _i ), (60) which concludes the proof. ∎ C.4 Dependent Rounding and Variance Reduction Throughout, Vart(X)≔t‖X−tX‖2Var_t(X) _t\|X-E_tX\|^2 denotes the scalar conditional variance of a vector. Let selection vector It∈0,1NI_t∈\0,1\^N indicate participating clients (Ii,t=(i∈St)I_i,t=I(i∈ S_t)). Algorithm 2 ensures: ∑i=1NIi,t=K,[Ii,t∣ℱt]=Kpi,Cov(Ii,t,Ij,t∣ℱt)≤0(∀i≠j). _i=1^NI_i,t=K, [I_i,t _t]=Kp_i, (I_i,t,I_j,t _t)≤ 0 (∀ i≠ j). (61) Lemma C.7 (Variance of Averaged Local Updates). Under dependent rounding with independent client execution, the aggregated update Mt=1K∑i=1NIi,tUi,tM_t= 1K _i=1^NI_i,tU_i,t satisfies: Var(Mt∣ℱt)≤1K∑i=1Npi[vi,t+2‖ai,t−a¯t‖2].Var(M_t _t)≤ 1K _i=1^Np_i [v_i,t+2\|a_i,t- a_t\|^2 ]. (62) Proof. Applying the law of total variance conditioned on filtration ℱtF_t: Vart(Mt)=t[Vart(Mt∣It)]+Vart(t[Mt∣It]).Var_t(M_t)=E_t [Var_t(M_t I_t) ]+Var_t (E_t[M_t I_t] ). (63) Conditional on selection ItI_t, client local trajectories are statistically independent: t[Vart(Mt∣It)]=t[1K2∑i=1NIi,tvi,t]=1K2∑i=1N(Kpi)vi,t=1K∑i=1Npivi,t.E_t [Var_t(M_t I_t) ]=E_t [ 1K^2 _i=1^NI_i,tv_i,t ]= 1K^2 _i=1^N(Kp_i)v_i,t= 1K _i=1^Np_iv_i,t. (64) For the second term, t[Mt∣It]=1K∑i=1NIi,tai,tE_t[M_t I_t]= 1K _i=1^NI_i,ta_i,t. Let C∈ℝN×NC ^N× N denote the covariance matrix of ItI_t conditional on ℱtF_t. Because ∑iIi,t=K _iI_i,t=K is strictly constant, ∑j=1NCij=0 _j=1^NC_ij=0, implying Cii=∑j≠i(−Cij)C_i= _j≠ i(-C_ij). Using negative correlation Cij≤0C_ij≤ 0 (i≠ji≠ j): Vart(1K∑i=1NIi,tai,t) _t ( 1K _i=1^NI_i,ta_i,t ) =1K2∑i=1N∑j=1NCij⟨ai,t,aj,t⟩ = 1K^2 _i=1^N _j=1^NC_ij a_i,t,a_j,t =12K2∑i≠j(−Cij)‖ai,t−aj,t‖2 = 12K^2 _i≠ j(-C_ij)\|a_i,t-a_j,t\|^2 ≤1K2∑i≠j(−Cij)(‖ai,t−a¯t‖2+‖aj,t−a¯t‖2) ≤ 1K^2 _i≠ j(-C_ij) (\|a_i,t- a_t\|^2+\|a_j,t- a_t\|^2 ) =2K2∑i=1NCii‖ai,t−a¯t‖2. = 2K^2 _i=1^NC_i\|a_i,t- a_t\|^2. (65) Since Ii,t∈0,1I_i,t∈\0,1\ is Bernoulli, Cii=Vart(Ii,t)=Kpi(1−Kpi)≤KpiC_i=Var_t(I_i,t)=Kp_i(1-Kp_i)≤ Kp_i. Substituting this into (65): Vart(1K∑i=1NIi,tai,t)≤2K∑i=1Npi‖ai,t−a¯t‖2.Var_t ( 1K _i=1^NI_i,ta_i,t )≤ 2K _i=1^Np_i\|a_i,t- a_t\|^2. (66) Combining (64) and (66) completes the proof. ∎ C.5 Multi-Step Trajectory Descent Guarantee The following proposition applies directly to the true E-step updates defined in (45), before bounding their bias and variance by model dimensions and heterogeneity. Proposition C.8 (Multi-Step Trajectory Descent Guarantee). Let Assumptions C.1–C.3 hold. Fix integers 1≤K≤N1≤ K≤ N and E,T≥1E,T≥ 1, and a constant local step size η>0η>0 with h=ηE≤1Lh=η E≤ 1L. Fix p before training, with ∑ipi=1 _ip_i=1 and 0≤pi≤1K0≤ p_i≤ 1K. At each round, draw St=DepRound(Kp)S_t=DepRound(Kp), independently of the conditionally independent client training streams, so that |St|=K|S_t|=K and Pr(i∈St∣ℱt)=Kpi (i∈ S_t _t)=Kp_i. Then, for Δ≔f(w0)−finf f(w_0)-f_ , 1T∑t=0T−1‖∇f(wt)‖2 1T _t=0^T-1E\|∇ f(w_t)\|^2 ≤2ΔhT+1T∑t=0T−1‖Bt(p)‖2 ≤ 2 hT+ 1T _t=0^T-1E\|B_t(p)\|^2 (67) +LhKT∑t=0T−1∑i=1Npi(vi,t+2∥ai,t−a¯t∥2), + LhKT _t=0^T-1E _i=1^Np_i (v_i,t+2\|a_i,t- a_t\|^2 ), where a¯t a_t and Bt(p)B_t(p) are defined in (46). Proof. By L-smoothness of global objective f(w)f(w) (Assumption C.1): f(wt+1)≤f(wt)+⟨∇f(wt),wt+1−wt⟩+L2‖wt+1−wt‖2.f(w_t+1)≤ f(w_t)+ ∇ f(w_t),w_t+1-w_t + L2\|w_t+1-w_t\|^2. (68) Substituting wt+1−wt=−hMtw_t+1-w_t=-hM_t and taking conditional expectations t[⋅]≔[⋅∣ℱt]E_t[·] [· _t]: tf(wt+1)≤f(wt)−h⟨gt,tMt⟩+Lh22(‖tMt‖2+Vart(Mt)).E_tf(w_t+1)≤ f(w_t)-h g_t,E_tM_t + Lh^22 (\|E_tM_t\|^2+Var_t(M_t) ). (69) Recall that tMt=∑ipiai,t=gt+Bt(p)E_tM_t= _ip_ia_i,t=g_t+B_t(p). Using the algebraic identity −2⟨a,b⟩+‖b‖2=−‖a‖2+‖b−a‖2-2 a,b +\|b\|^2=-\|a\|^2+\|b-a\|^2: −h⟨gt,tMt⟩+h2‖tMt‖2=−h2‖gt‖2+h2‖Bt(p)‖2.-h g_t,E_tM_t + h2\|E_tM_t\|^2=- h2\|g_t\|^2+ h2\|B_t(p)\|^2. (70) Since effective step size satisfies Lh≤1Lh≤ 1, we have Lh22≤h2 Lh^22≤ h2. Substituting (70) into (69): tf(wt+1)≤f(wt)−h2‖gt‖2+h2‖Bt(p)‖2+Lh22Vart(Mt).E_tf(w_t+1)≤ f(w_t)- h2\|g_t\|^2+ h2\|B_t(p)\|^2+ Lh^22Var_t(M_t). (71) Rearranging and bounding Vart(Mt)Var_t(M_t) via Lemma C.7: ‖gt‖2≤2(f(wt)−tf(wt+1))h+‖Bt(p)‖2+LhK∑i=1Npi(vi,t+2‖ai,t−a¯t‖2).\|g_t\|^2≤ 2(f(w_t)-E_tf(w_t+1))h+\|B_t(p)\|^2+ LhK _i=1^Np_i (v_i,t+2\|a_i,t- a_t\|^2 ). (72) Taking total expectations, summing over t=0,…,T−1t=0,…,T-1, telescoping the objective values ∑t=0T−1(f(wt)−f(wt+1))=f(w0)−f(wT)≤f(w0)−finf=Δ _t=0^T-1(Ef(w_t)-Ef(w_t+1))=f(w_0)-Ef(w_T)≤ f(w_0)-f_ = , and dividing by T completes the proof of Proposition C.8. ∎ C.6 Formal Non-Convex Convergence Guarantee The following theorem gives the formal version of Theorem 4.1, including the sampling, horizon, and parameter requirements. Let aH,ζ≥0a_H,ζ≥ 0 be the constants in Assumption C.4, and define Mp≔∑i=1NpidiZO(diZO+2)(diZO+4),χ2(p)≔N∑i=1N(pi−1N)2.M_p _i=1^Np_id_i^ZO(d_i^ZO+2)(d_i^ZO+4), χ^2(p) N _i=1^N (p_i- 1N )^2. (73) Thus MpM_p controls the finite-difference remainder, and χ2(p)χ^2(p) measures departure from uniform participation. Theorem C.9 (Non-Convex Convergence of HO-FL with Fixed Sampling). Suppose Assumptions C.1–C.4 hold, with σi2≤σ2 _i^2≤σ^2 and Lis≤LsL_i^s≤ L_s for all clients. Fix integers 1≤K≤N1≤ K≤ N and E,q,T≥1E,q,T≥ 1, client boundaries, and a probability vector p satisfying ∑ipi=1 _ip_i=1 and ϵN≤pi≤1K εN≤ p_i≤ 1K for a fixed ϵ∈(0,1]ε∈(0,1]. Use the local updates (7) and the dependent-rounding and independence conditions of Proposition C.8. The following explicit constants are sufficient: C0 C_0 ≔1+(30+12ϵ)[(32N+8σ2)L2+2Ls2]+(54N+12σ2)L, 1+ (30+ 12ε ) [(32N+8σ^2)L^2+2L_s^2 ]+(54N+12σ^2)L, (74) C⋆ C_ ≔max1, 8L2(1+2Nϵ),(2NLϵ)2, 144C02(1+aH)2. \1,\;8L^2 (1+ 2Nε ),\; ( 2NLε )^2,\;144C_0^2(1+a_H)^2 \. Assume the sampling skew and training horizon satisfy aHχ2(p)≤112,T≥C⋆(1+1q∑i=1NpidiZO).a_Hχ^2(p)≤ 112, T≥ C_ (1+ 1q _i=1^Np_id_i^ZO ). (75) Choose the common local step size and perturbation radius as η η =1ET(1+1q∑i=1NpidiZO), = 1E T (1+ 1q _i=1^Np_id_i^ZO ), (76) 0<μ 0<μ ≤1Mp(1+1q∑i=1NpidiZOT)14if Mp>0. ≤ 1 M_p ( 1+ 1q _i=1^Np_id_i^ZOT ) 14 M_p>0. All clients use the same q and μ; if Mp=0M_p=0, all clients use FO and no perturbation radius is needed. Then, with Δ=f(w0)−finf =f(w_0)-f_ and Crate≔4Δ+(8+6ζ2)C0C_rate 4 +(8+6ζ^2)C_0, 1T∑t=0T−1‖∇f(wt)‖2≤Crate1+1q∑i=1NpidiZOT+6ζ2N∑i=1N(pi−1N)2. 1T _t=0^T-1E\|∇ f(w_t)\|^2≤ C_rate 1+ 1q _i=1^Np_id_i^ZOT+6ζ^2N _i=1^N (p_i- 1N )^2. (77) In particular, this yields the two-term O rate in (11). Its implicit constants can depend on N,ϵ,L,Ls,σ,aH,ζN,ε,L,L_s,σ,a_H,ζ and Δ , but not on E,T,K,q,μ,diZOE,T,K,q,μ,d_i^ZO or pip_i. Proof. We first verify that the stated horizon and parameter choices imply every stability condition used below. Let κmax=maxiκi _ = _i _i and τ=η2(E−1)2τ=η^2(E-1)^2. Coverage and diZO+1≤2diZOd_i^ZO+1≤ 2d_i^ZO for a ZO client imply κmax≤2Nϵ(1+1q∑ipidiZO). _ ≤ 2Nε (1+ 1q _ip_id_i^ZO ). Substituting (76) and using (75) gives LηE Lη E ≤LC⋆≤1, ≤ L C_ ≤ 1, (78) 4(1+κmax)L2τ 4(1+ _ )L^2τ ≤4L2(1+2Nϵ)C⋆≤12, ≤ 4L^2(1+ 2Nε)C_ ≤ 12, LηκmaxK Lη _ K ≤2NLϵC⋆≤1. ≤ 2NLε C_ ≤ 1. Thus Proposition C.8 and Lemma C.6 apply. Define the proof coefficient Ap(η)≔η(E+1+1q∑ipidiZO)K+η2(E−1)2(1+1q∑ipidiZO).A_p(η) η (E+1+ 1q _ip_id_i^ZO )K+η^2(E-1)^2 (1+ 1q _ip_id_i^ZO ). (79) Since E,K≥1E,K≥ 1, the chosen step size bounds the local-drift term by 1T 1T and gives Ap(η) A_p(η) ≤21+1q∑ipidiZOT+1T ≤ 2 1+ 1q _ip_id_i^ZOT+ 1T (80) ≤31+1q∑ipidiZOT≤3C⋆. ≤ 3 1+ 1q _ip_id_i^ZOT≤ 3 C_ . By (74), this ensures C0(1+aH)Ap(η)≤14C_0(1+a_H)A_p(η)≤ 14. Next, fix a round and suppress t. Write g=∇f(wt)g=∇ f(w_t), gi=∇fi(wt)g_i=∇ f_i(w_t), H=H(wt)H=H(w_t), G2=‖g‖2+HG_2=\|g\|^2+H, and V=Var(Mt∣ℱt)V=Var(M_t _t). Using Qi,DiQ_i,D_i from Lemma C.6, define Qp=∑ipiQiQ_p= _ip_iQ_i, D¯=∑ipiDi D= _ip_iD_i, νp=∑ipiνi _p= _ip_i _i and ρp2=∑ipiρi2 _p^2= _ip_i _i^2. Since ‖gi−g‖2≤NH\|g_i-g\|^2≤ NH and ∑ipi(1+κi)≤2(1+1q∑ipidiZO) _ip_i(1+ _i)≤ 2(1+ 1q _ip_id_i^ZO), the estimator moments imply Qp Q_p ≤4N(1+1q∑ipidiZO)G2, ≤ 4N (1+ 1q _ip_id_i^ZO )G_2, νp _p ≤4σ2(1+1q∑ipidiZO)+Ls22μ2Mp,ρp2≤Ls2μ2Mp, ≤ 4σ^2 (1+ 1q _ip_id_i^ZO )+ L_s^22μ^2M_p, _p^2≤ L_s^2μ^2M_p, D¯ D =2τ(4Qp+νp). =2τ(4Q_p+ _p). (81) Cauchy–Schwarz gives ‖∑i(pi−1N)gi‖2≤χ2(p)H. \| _i (p_i- 1N )g_i \|^2≤χ^2(p)H. (82) Combining this with the bias decomposition and Lemma C.6 yields ‖B(p)‖2≤3χ2(p)H+3L2D¯+3ρp2.\|B(p)\|^2≤ 3χ^2(p)H+3L^2 D+3 _p^2. (83) For the participation variance, let ri=ai−gir_i=a_i-g_i and r¯=1N∑iri r= 1N _ir_i. The coverage condition gives ‖r¯‖2≤1ϵ∑ipi‖ri‖2\| r\|^2≤ 1ε _ip_i\|r_i\|^2, so ∑ipi‖ai−a¯‖2≤3NH+6(1+1ϵ)(L2D¯+ρp2). _ip_i\|a_i- a\|^2≤ 3NH+6 (1+ 1ε )(L^2 D+ _p^2). (84) Together with (52), Lemma C.7, and the verified stability bounds (78), this gives ‖B(p)‖2+LhV \|B(p)\|^2+LhV ≤3χ2(p)H+(30+12ϵ)L2D¯+(18+12ϵ)ρp2 ≤ 3χ^2(p)H+ (30+ 12ε )L^2 D+ (18+ 12ε ) _p^2 (85) +6NLηEKH+12LηKQp+3LηKνp. + 6NLη EKH+ 12LηKQ_p+ 3LηK _p. Substitute (81). Since (78) also implies L2τ≤18L^2τ≤ 18 and LηK≤1 LηK≤ 1, the constant C0C_0 in (74) bounds all resulting coefficients, yielding ‖B(p)‖2+LhV≤3χ2(p)H+C0Ap(η)(‖g‖2+H+1)+C0μ2Mp.\|B(p)\|^2+LhV≤ 3χ^2(p)H+C_0A_p(η)(\|g\|^2+H+1)+C_0μ^2M_p. (86) For all-FO training, the terms involving μ2Mpμ^2M_p are understood as zero. By (75), 3aHχ2(p)≤143a_Hχ^2(p)≤ 14, and the first part of the proof already established C0(1+aH)Ap(η)≤14C_0(1+a_H)A_p(η)≤ 14. Substituting the affine heterogeneity bound into (86) therefore contributes at most 12‖g‖2 12\|g\|^2 to the right-hand side of Proposition C.8. Absorbing this term gives 1T∑t=0T−1‖gt‖2≤4ΔηET+2C0(1+ζ2)Ap(η)+2C0μ2Mp+6ζ2χ2(p). 1T _t=0^T-1E\|g_t\|^2≤ 4 η ET+2C_0(1+ζ^2)A_p(η)+2C_0μ^2M_p+6ζ^2χ^2(p). (87) Finally, 1ηET=(1+1q∑ipidiZO)/T 1η ET= (1+ 1q _ip_id_i^ZO)/T, while (76) bounds μ2Mpμ^2M_p by the same quantity. Substitution into (87), together with (80), gives (77) with Crate=4Δ+(8+6ζ2)C0C_rate=4 +(8+6ζ^2)C_0. ∎ The step size controls local drift directly, yielding a rate that isolates dimension and sampling effects without claiming speedup in E or K. Uniform participation removes the sampling-bias term; non-uniform participation trades a smaller sampling-weighted ZO dimension against residual data heterogeneity. Appendix D Dimension-Aware Sampling: Derivations and Practical Guidance This appendix supplies the derivations behind Section 5 and a practical procedure for choosing the fixed sampling preference β. We first derive the scalar relaxation and its constrained solution, then quantify the role of β and describe how to select it before training. D.1 Derivation of the Sampling Relaxation The statistical reference problem. Condition on the history before a round, and treat ai,vi,ga_i,v_i,g as fixed conditional statistics. Set =[a1,…,aN]A=[a_1,…,a_N] and ci=LhK[vi+2‖ai−a¯‖2]c_i= LhK[v_i+2\|a_i- a\|^2]. The sampling-dependent part of Proposition C.8 is J(p)=‖p−g‖2+c⊤p=p⊤⊤p−2g⊤p+‖g‖2+c⊤p.J(p)=\|Ap-g\|^2+c p=p A Ap-2g Ap+\|g\|^2+c p. (88) Its Hessian is 2⊤⪰02A A 0, establishing convexity of (13). This is an optimum of the current conditional bound; the practical policy uses a single probability vector fixed before training. Scalar relaxation. Let ui=1Nu_i= 1N and Ha=1N∑i‖ai−a¯‖2H_a= 1N _i\|a_i- a\|^2. Since ∑i(pi−ui)=0 _i(p_i-u_i)=0, p−g -g =∑i(pi−ui)(ai−a¯)+(a¯−g), = _i(p_i-u_i)(a_i- a)+( a-g), ‖p−g‖2 \|Ap-g\|^2 ≤2‖∑i(pi−ui)(ai−a¯)‖2+2‖a¯−g‖2 ≤ 2 \| _i(p_i-u_i)(a_i- a) \|^2+2\| a-g\|^2 ≤2‖p−u‖2∑i‖ai−a¯‖2+2‖a¯−g‖2. ≤ 2\|p-u\|^2 _i\|a_i- a\|^2+2\| a-g\|^2. (89) Thus J(p)≤A‖p−u‖2+c⊤p+2‖a¯−g‖2J(p)≤ A\|p-u\|^2+c p+2\| a-g\|^2, with A=2NHaA=2NH_a. The last term is independent of p, yielding (14). This removes g from the sampling optimization, while its effect remains in the convergence analysis through a¯−g a-g. The relaxation preserves client-level noise and update dispersion but discards directional cancellation and the cross term with a¯−g a-g. KKT conditions and the clipped solution. For A>0A>0, first consider K<NK<N and ϵ<1ε<1, and write p¯=ϵN p= εN and p¯=1K p= 1K. The Lagrangian is (p,λ,α,ω)=A‖p−u‖2+c⊤p+λ(∑ipi−1)+∑iαi(p¯−pi)+∑iωi(pi−p¯),J(p,λ,α,ω)=A\|p-u\|^2+c p+λ ( _ip_i-1 )+ _i _i( p-p_i)+ _i _i(p_i- p), where αi,ωi≥0 _i, _i≥ 0. Stationarity and complementary slackness require 2A(pi−ui)+ci+λ−αi+ωi 2A(p_i-u_i)+c_i+λ- _i+ _i =0, =0, (90) αi(pi−p¯)=0,ωi(p¯−pi) _i(p_i- p)=0, _i( p-p_i) =0. =0. At an interior coordinate, both inequality multipliers vanish, giving pi=ui−ci+λ2Ap_i=u_i- c_i+λ2A. At the lower bound, stationarity implies that this unconstrained value is at most p¯ p; at the upper bound it is at least p¯ p. Consequently, pi⋆=clip[ϵ/N, 1/K](1N−ci+λ2A).p_i =clip_[ε/N,\,1/K] ( 1N- c_i+λ2A ). (91) Strict convexity makes p⋆p unique. The sum of its coordinates is continuous and nonincreasing in λ. A valid bisection bracket is λlo=mini2A(1N−1K)−ci,λhi=maxi2A(1N−ϵN)−ci. _lo= _i \2A ( 1N- 1K )-c_i \, _hi= _i \2A ( 1N- εN )-c_i \. At these endpoints, the coordinate sums are NK≥1 NK≥ 1 and ϵ≤1ε≤ 1, respectively. When A=0A=0, the objective is linear: initialize every coordinate at ϵN εN and allocate the remaining mass in increasing order of cic_i, up to 1K 1K; tied costs can share the allocation. If K=NK=N or ϵ=1ε=1, the feasible set is simply u\u\. D.2 Choosing β in Practice HO-FL uses β as a single preference parameter, selected before the training run and held fixed. The following guidelines combine scenario knowledge, a small validation budget, and desired participation coverage. Starting from scenario knowledge. When client data are broadly representative of the same population and ZO dimensions differ substantially, β=0.75β=0.75 is a useful starting point. For intermediate heterogeneity, or limited prior knowledge, start from β=0.5β=0.5. With pronounced client-specific data skew, start from β=0.25β=0.25; β=0β=0 remains the uniform reference. Resource–data correlation matters as well: if the preferred clients represent only a narrow part of the population, a smaller preference is appropriate even when their dimension advantage is large. These values are empirical starting points, rather than universal heterogeneity thresholds. Table 4 illustrates the pattern on Qwen/AG News: 0.750.75 performs best under IID and the two milder skews, whereas 0.250.25 performs best at α=0.1α=0.1. Further IID comparisons appear in Appendix E.4. A small-budget validation procedure. Without reliable scenario knowledge, we recommend comparing β∈0,0.25,0.5,0.75β∈\0,0.25,0.5,0.75\ in short pilot runs on a representative subset of the target task. Preserve the client partition structure, resource tiers, model boundaries, and any resource–data correlation. Evaluate on held-out validation data spanning all resource tiers, using the mean client validation loss to match the uniformly weighted client objective. Select the best candidate, favoring a smaller β when validation differences are inconclusive, then fix it for the formal run. Federated validation and such successive-halving procedures are established approaches to hyperparameter selection (Khodak et al., 2021). This is a recommended deployment procedure; the experimental configurations and sensitivity sweeps used in this paper are reported in Appendix E. Choosing a coverage budget. A complementary criterion requires no gradient or heterogeneity estimates. For a nonpreferred client i∉Ri∉ R, exact inclusion probabilities imply [∑t=0T−1i∈St]=TKpi=KTN(1−β).E\! [ _t=0^T-11\i∈ S_t\ ]=TKp_i= KTN(1-β). (92) If the desired minimum expected participation count is m∈[0,KT/N]m∈[0,KT/N], restrict candidates to β≤1−mNKT.β≤ 1- mNKT. (93) For our N=30N=30, K=6K=6, and T=160T=160 setting, β=0.25,0.5,0.75β=0.25,0.5,0.75 correspond to 24,16,824,16,8 expected participations per nonpreferred client. This is an expectation, rather than a guaranteed count in every run. It translates the preference parameter into a concrete participation budget: choose an acceptable coverage level, then use scenario knowledge or validation to select among the remaining candidates. Once selected, β and the resource ranking remain fixed, so this procedure adds no online moment estimation to the training algorithm. Appendix E Experimental Details E.1 Datasets and Benchmark Configurations We evaluate all frameworks across five representative language understanding and generation benchmarks. Table 6 summarizes their task types, sample sizes, and evaluation metrics: • SST-2 (Socher et al., 2013): The Stanford Sentiment Treebank binary classification task. Evaluated on the full official validation set (872 examples). • BoolQ (Clark et al., 2019): A reading comprehension benchmark of natural yes/no questions paired with Wikipedia passages. Evaluated on the full official validation set (3,270 examples). • SciQ (Welbl et al., 2017): A crowdsourced multiple-choice science examination benchmark spanning physics, chemistry, and biology. Evaluated on the official test set (1,000 examples). • AG News (Zhang et al., 2015): A 4-class topic classification dataset (World, Sports, Business, Sci/Tech). We evaluate IID and non-IID settings with 4,000 training examples per client to control for quantity skew (Li et al., 2022). Non-IID class proportions follow Dir(α4)Dir( 1_4), α∈1,0.5,0.1α∈\1,0.5,0.1\; exhausted class pools are replenished by sampling with replacement, allowing duplicates within and across clients. Evaluation uses the standard test set (7,600 examples). • SQuAD 1.1 (Rajpurkar et al., 2016): A challenging token-level extractive question answering benchmark. We adopt greedy auto-regressive decoding with a maximum prompt length of 512 tokens and a generation cap of 64 tokens. We report official Exact Match (EM) and Macro-F1 metrics on all 10,570 validation examples. Table 6: Benchmark dataset statistics and evaluation protocols. Dataset Task Train Eval Partition Metric SST-2 Sentiment Classification 67,349 872 IID Accuracy (%) BoolQ Reading Comprehension (Y/N) 9,427 3,270 IID Accuracy (%) SciQ Science QA (4 choices) 11,679 1,000 IID Accuracy (%) AG News Topic Classification (4 classes) 120,000 7,600 IID / Dir(α)Dir(α) Accuracy (%) SQuAD 1.1 Extractive QA 87,599 10,570 IID EM / F1 (%) E.2 Model Architectures and LoRA Configuration All pretrained foundation models are instantiated using their official HuggingFace checkpoints: facebook/opt-125m (12 transformer layers, hidden dimension 768), Qwen/Qwen2.5-1.5B (28 layers, hidden dimension 1536), and HuggingFaceTB/SmolLM3-3B (36 layers, hidden dimension 2048). The backbone network weights are fully frozen and cast to BF16 precision. Low-Rank Adaptation (LoRA) (Hu et al., 2022) is injected exclusively into the query and value projection matrices (Wq,WvW_q,W_v) across all attention modules: • LoRA Hyperparameters: We configure adapter rank r=8r=8, LoRA scaling factor αlora=16 _lora=16, and zero dropout across all experiments. LoRA matrices A and B are stored in FP32 precision to guarantee numerical stability during zeroth-order gradient perturbations and finite-difference evaluations. • Optimizer Configuration: Local training uses AdamW (Loshchilov and Hutter, 2019) with decoupled weight decay 5×10−45× 10^-4, learning rate 10−510^-5, momentum parameters (β1,β2)=(0.9,0.999)( _1, _2)=(0.9,0.999), and stabilizer ϵ=10−8ε=10^-8. Optimizer momentum buffers are reset upon each client participation. E.3 Baseline Implementation Details We compare HO-FL against the following state-of-the-art baselines using matched data partitions and seeds at the same number of communication rounds: • FedAvg (McMahan et al., 2017): Executes standard backpropagation (FO) on all trainable LoRA adapters across the entire model. All selected clients update all adapters and transmit their parameter deltas. • FedProx (Li et al., 2020): Adds a proximal regularization term μprox2‖w−wt‖2 _prox2\|w-w_t\|^2 to the local objective to combat client drift under non-IID data distributions. We set μprox=10−3 _prox=10^-3. • DeComFL (Li et al., 2025): A state-of-the-art pure zeroth-order federated optimization method. All trainable parameters are perturbed via randomized Gaussian directions (q=2,μ=10−3q=2,μ=10^-3). Communication cost is dimension-free via shared seed generation. We adapt DeComFL to support local AdamW updates with deterministic replay. • HO-FL (β=0β=0): Executes hybrid-order local updates identical to HO-FL across the heterogeneous client tiers, but samples participating clients uniformly at random (pi=1/Np_i=1/N, equivalent to setting β=0β=0). • HO-FL (Proposed): Employs the dimension-aware client sampling policy p=(1−β)u+βrp=(1-β)u+β r with dependent rounding St=DepRound(Kp)S_t=DepRound(Kp), using β=0.75β=0.75 for IID benchmarks and exploring β∈[0,0.75]β∈[0,0.75] on Dirichlet label partitions. Training and evaluation. All models use AdamW with learning rate 10−510^-5, moment coefficients (0.9,0.999)(0.9,0.999), numerical stabilizer 10−810^-8, and weight decay 5×10−45× 10^-4; optimizer states reset at each client participation. FedProx uses a proximal coefficient of 10−310^-3. The frozen backbone uses BF16, and the query/value LoRA adapters use FP32 with rank 8, scaling parameter 16, and zero dropout. The default configuration has N=30N=30, K=6K=6, batch size 16, E=5E=5, q=2q=2, μ=10−3μ=10^-3, and 160 rounds. SST-2 and BoolQ use their full validation sets of 872 and 3,270 examples. SciQ uses the official 1,000-example test set. SQuAD uses greedy generation with a prompt cap of 512 tokens and an output cap of 64 tokens; Table 3 reports official v1.1 EM/F1 on all 10,570 validation examples. Intermediate SQuAD curves use a fixed 1,024-example subset. We report the final model throughout. Seeds and sampling. All entries in Table 2, the Qwen/AG News sensitivity study, and the local-training ablations use three seeds, reporting means and sample SDs. All SQuAD results use these three seeds. All SQuAD sampling comparisons use three matched seeds. Within each seed, comparisons preserve the client data, initialization, and resource assignment. The preferred set R is fixed using the smallest actual trainable ZO dimensions divided by q, with an independent seeded tie-breaking permutation. HO-FL applies dependent rounding to KpKp and directly averages the selected models. E.4 IID sampling preference Table 7 compares sampling strengths using seed 42; Table 8 uses three seeds. Each comparison preserves the configuration within each seed. We use β=0.75β=0.75 for the IID main results. Relative to 0.50.5, it preserves SST-2 accuracy for both OPT and Qwen, preserves OPT/BoolQ accuracy, and improves Qwen/BoolQ and Qwen/SciQ by 0.55 and 0.30 points. OPT/SciQ decreases by 0.20 points while remaining near chance. On SQuAD, 0.750.75 improves both EM and F1 over 0.50.5 for all three models in every seed; the mean F1 gains are 1.27, 0.30, and 0.27 points for OPT, Qwen, and SmolLM3. The three-seed Qwen/AG News IID comparison in Table 4 also favors 0.750.75 over 0.50.5 in every seed. Table 7: IID sampling sensitivity. Final accuracy (%). All columns use HO-FL. Model Task β=0β=0 0.250.25 0.50.5 0.750.75 OPT SST-2 86.01 86.47 86.47 86.47 BoolQ 57.98 57.92 58.84 58.84 SciQ 24.20 23.90 24.50 24.30 Qwen SST-2 91.28 91.86 92.78 92.78 BoolQ 76.54 77.74 78.29 78.84 SciQ 89.90 89.80 90.00 90.30 Table 8: SQuAD sampling sensitivity. Final full-validation EM/F1 (%), reported as mean(SD)mean_(SD) accross three seeds. Model Metric β=0β=0 0.50.5 0.750.75 OPT EM 37.16(0.26)37.16_(0.26) 40.68(0.23)40.68_(0.23) 41.87(0.32)41.87_(0.32) F1 47.24(0.21)47.24_(0.21) 51.16(0.19)51.16_(0.19) 52.43(0.20)52.43_(0.20) Qwen EM 76.67(0.05)76.67_(0.05) 78.32(0.18)78.32_(0.18) 78.71(0.19)78.71_(0.19) F1 85.34(0.03)85.34_(0.03) 86.61(0.11)86.61_(0.11) 86.92(0.10)86.92_(0.10) Smol EM 77.80(0.11)77.80_(0.11) 79.23(0.06)79.23_(0.06) 79.54(0.20)79.54_(0.20) F1 87.25(0.04)87.25_(0.04) 88.30(0.06)88.30_(0.06) 88.57(0.13)88.57_(0.13) Stronger preference increases the participation of clients with larger FO-trained segments. With three equally sized resource tiers, their expected per-round counts are 4/1/14/1/1 at β=0.5β=0.5 and 5/0.5/0.55/0.5/0.5 at 0.750.75, ordered from the smallest to the largest ZO segment. The choice of 0.750.75 therefore favors task performance while retaining the per-device memory benefit of hybrid training. The heterogeneity study reports the full β grid, and the local-training ablations below keep their fixed β=0.5β=0.5 setting. E.5 Memory profiles and local participation simulation Profiles use Qwen2.5-1.5B, batch size 16, E=5E=5, and q=2q=2, with FP32 query/value LoRA adapters and a frozen BF16 backbone. We use the maximum peak allocated bytes from three isolated measurements of a complete local participation, including optimizer state initialization. Only sequence length 1024 is displayed; FedProx is omitted from the memory comparison. The right panel of Figure 2 is generated locally, without training a model. It retains the seed-42 client cut assignment and preferred set from the Qwen experiments, with ten clients each at cuts 8/16/24. We set pi=(1−β)/30+βi∈R/6p_i=(1-β)/30+ 1\i∈ R\/6 and simulate 160 rounds of fixed-size dependent rounding for each β∈0,0.25,0.5,0.75β∈\0,0.25,0.5,0.75\. Each round uses a fresh random processing order and randomness derived from the simulation seed and round index. Client profiles are summed for the six selected participants and averaged over rounds. The corresponding expectation is 6∑ipimi(L)6 _ip_im_i(L); finite-round averages need not equal this expectation exactly. FedAvg and DeComFL use homogeneous profiles, so their round totals are six times their single-client requirement. E.6 Local training ablations We vary one parameter at a time on OPT-125M/SST-2 IID with HO-FL, fixing β=0.5β=0.5 and using three random seeds. The default E=5E=5, q=2q=2, μ=10−3μ=10^-3 configuration is shared across the sweeps. Table 9 reports final full-validation accuracy and task loss. Table 9: Local-training sensitivity on OPT/SST-2 IID. HO-FL with β=0.5β=0.5; mean ± sample SD across three seeds. All configurations process 76,800 training examples. The default is shared across the three sweeps. Configuration E q μ Rounds Accuracy (%) Validation loss Default 5 2 10−310^-3 160 86.39±0.0786.39± 0.07 0.3498±0.00380.3498± 0.0038 E=1E=1 1 2 10−310^-3 800 86.47±0.4186.47± 0.41 0.3489±0.00330.3489± 0.0033 E=10E=10 10 2 10−310^-3 80 86.58±0.3086.58± 0.30 0.3406±0.00080.3406± 0.0008 q=1q=1 5 1 10−310^-3 160 86.51±0.1886.51± 0.18 0.3498±0.00360.3498± 0.0036 q=4q=4 5 4 10−310^-3 160 86.47±0.0086.47± 0.00 0.3498±0.00380.3498± 0.0038 q=8q=8 5 8 10−310^-3 160 86.47±0.1186.47± 0.11 0.3496±0.00410.3496± 0.0041 μ=10−4μ=10^-4 5 2 10−410^-4 160 86.31±0.0786.31± 0.07 0.3496±0.00360.3496± 0.0036 μ=10−2μ=10^-2 5 2 10−210^-2 160 86.35±0.0086.35± 0.00 0.3493±0.00400.3493± 0.0040 Local steps E. For E∈1,5,10E∈\1,5,10\, we use 800, 160, and 80 rounds, respectively, so each run processes 76,800 training examples. Mean accuracy is 86.47%, 86.39%, and 86.58%. In this range, increasing local work preserves performance while reducing the number of communication rounds; E=10E=10 also yields a lower final validation loss. Direction count q. For q∈1,2,4,8q∈\1,2,4,8\, mean accuracy is 86.51%, 86.39%, 86.47%, and 86.47%. Increasing the direction count does not produce a monotonic performance gain in this setting. A small number of directions is sufficient, supporting the default q=2q=2 with its lower local computation cost. Perturbation radius μ. The radii 10−410^-4, 10−310^-3, and 10−210^-2 yield 86.31%, 86.39%, and 86.35% mean accuracy. The small variation across this range supports retaining μ=10−3μ=10^-3. Together, these sweeps show that performance is stable across the tested local-training parameters.