Nonlinear-Gain Distributed Zeroth-Order Optimization for Networked Black-Box Control

This letter studies distributed stochastic optimization over a peer-to-peer network when agents can query only zeroth-order function values, proposing ZOOM-PB, a coordinate-sampling method that blends each local ZO estimate with a fractional-power response while maintaining only a primal state. The raw estimate is retained as a linear anchor, and the nonlinear mixing weight is coupled to the optimization stepsize. This design is motivated by a basic obstruction: transforming heterogeneous or noisy local estimates before averaging can reverse the network direction. The paper bounds that nonlinear residual directly from the raw oracle assumptions instead of imposing an aggregate-alignment condition. With a smooth stochastic-function oracle and a connected graph, ZOOM-PB attains the nonconvex stationarity order O(p/(nT)) and a Polyak–Łojasiewicz statistical term of order O(p/(nT)), after an explicit initialization transient. Numerical examples compare ZOOM-PB with seven distributed ZO baselines under matched query and message budgets.

Distributed optimization is a basic tool for networked control, machine learning, robotic swarms, and sensor networks, where agents cooperatively solve

minp f (x) ≜

x∈R

1X

fi (x),

n i=1

fi (x) = Eξi [Fi (x, ξi)].

Many first-order distributed methods assume that each agent can compute or estimate ∇fi (x) directly. This assumption is restrictive in black-box control and simulation-based optimization, where agents may observe only noisy function values, such as concentration readings in source seeking [1] or loss values in black-box learning [2]–[5]. Zeroth-order (ZO) methods replace gradients by finite-difference estimates and have been extensively studied in centralized stochastic optimization [6]–[9]. Distributed ZO methods further couple gradient estimation with network consensus [10]–[16]. However, most existing distributed stochastic ZO algorithms follow one of two paths: they either

Pγ,τ (g) = τ 1−γ sgn(g)gγ,

(2)

applied componentwise. Its gain ratio is Pγ,τ (g)

= g

τ

g

1−γ

Aγ,τ,β (g) = (1 − β)g + βPγ,τ (g),

(3)

β ∈ [0, 1].

This work was partially supported by the National Natural Science Foundation of China under grant 52307120 and the National Social Science Fund of China under grant 24CTJ010.

The map is evaluated locally; β = 0 gives the raw ZO direction and β = 1 gives the pure powerball direction.

ZOOM-PB combines sampled coordinate queries with a one-state recursion: each agent maintains and broadcasts only xi,k. Local oracle calls run in parallel, but each agent still queries only its own objective. Second, we show that the nonlinear residual in (4) is a controlled perturbation under the same raw stochastic-oracle assumptions used for the linear direction. Choosing βk2 = O(ηk /n) preserves the stated nonconvex and PL statistical orders without assuming that the averaged purepowerball direction is aligned. The proof also bounds the nonlinear forcing in the consensus recursion. Third, matchedbudget experiments separate gain shaping from estimator and state choices across coordinate, spherical, and Gaussian distributed-ZO baselines.

Proposition 1. For γ ∈ (0, 1), τ > 0, β ∈ (0, 1], and g ≠ 0,

Aγ,τ,β preserves the sign of g, and

Aγ,τβ (g)

= (1 − β) + β

g

It amplifies components below τ, leaves g = τ unchanged, and attenuates components above τ. Moreover, A1,τβ (g) = g. Proof. Both terms in (4) have the sign of g. Dividing their convex combination by g gives the displayed ratio.

Algorithm 1 ZOOM-PB

1: Input: α > 0, γ ∈ [1/2, 1], τ > 0, steps ηk, weights βk ∈ [0, 1], radii δi,k.

2: Initialize: xi,0 ∈ Rp for all i ∈ [n].

3: for k = 0, 1,... do

4: for each agent i in parallel do

5: Exchange xi,k with neighbors j ∈ Ni.

6. Form gi,k by (5) or (6).

7. Set si,k = Aγ,τβk gi,k

8. Update xi,k+1 = xi,k − αnX Lij xj,k − ηkk si,k.

9: end for

10: end for

Remark 2. For γ = 1 or βk = 0, ri,k = 0 and the analysis reduces to the usual ZOOM proof. For γ < 1, schedule (11) controls the possibly misaligned nonlinear residual. The pair (γ, τ) sets the nonlinear response. Smaller γ increases both weak-component actuation and sensitivity to noise near zero; τ records the normalization at which amplification changes to attenuation. The weight βk controls departure from the raw ZO direction. Equation (11) gives the admissible range βk ∈ [0, β̄k]. The condition 16Rg βd2 ≤ 1 makes the worst-case residual subordinate to raw descent; the stepsize-dependent cap preserves the rate order. None of these parameters changes the query or message budget, and βk requires no stored vector state. They should be reported with the objective scaling and selected by a validation or noise-sensitivity study. The radius schedules in the theorems correspond to the common-random-number oracle; independent measurement noise requires the modification in Remark 1.

Theorem 1. Suppose Assumptions 1–3 hold, p/nc ≤ κc, W0 ≤ CW n, and

2, λn (L) r

n, ηk = η = pT

0<α<

0 ≤ βk ≤ β̄k, δi,k

E[χT] ≤ CρT χ0 +

Cp

CW0 ta0

+. (T + t0)2

n(T + t0)a+2

Proof. The PL inequality and (12) imply, for a constant cV > 0,

E[Wk+1] ≤ (1 − cV ηk)E[Wk] + C0 pηk squared + Cδ np squared ηk δ i,k squared. The chosen radius makes the last term at most Cpη i squared; the nonlinear residual has the same order by (11). Thus

aCp

uk+1 ≤ 1 −,uk = E[Wk].

The product estimate in Appendix D gives

a t0 Cp,uT ≤ C u0 + T + t0

T + t0

κδ <= 3/4 1/4. p n (k + 1)1/4

If T ≥ n3 /(pη̄ squared), then nη ≤ η̄, the restriction used in Appendix B.

Theorem 2. Suppose Assumptions 1–4 hold and p/nc ≤ κc. Choose

r

κη ηk =,0 ≤ βk ≤ β̄k, δi,k ≤ κδ, k + t0 np which proves (21). Finally, (17),∥∇f (x)∥ ≤ 2Lf (f (x − f ∗)), and this finite-time bound give a geometrically stable recursion whose forcing is Cp/(k + t0)2 + CW0 ta0 /[n(k + t0)a+2]. Its convolution proves (22).

The consensus decrease in (12) dominates cηk χk for all sufficiently small ηk. Hence there is a cV > 0, depending on ν and the graph constants, such that

• Cδ np squared ηk δ i,k squared.

E[χT] ≤ C. T

T

(42)

Finally, T ≥ n /(pη̄), the restriction used in Appendix B.

The paper concludes that ZOOM-PB combines coordinate function queries, a componentwise nonlinear gain, and a one-state recursion. Retaining the raw ZO direction as a linear anchor avoids an aggregate-alignment assumption: the possibly misaligned nonlinear component is controlled as a vanishing perturbation. The resulting method has the stated nonconvex order and PL statistical term without an auxiliary tracking state. The broad classification comparison gives comparable endpoints. In the weak-signal UAV sweep, ZOOM-PB uses 43–74% fewer queries at the two held-out weakest scales than ZOOM, ZOD-PA, and ZOD-PDA under fixed per-method tuning; the noise sweep shows why this is not uniform superiority. Extending the theory to independent measurement noise and adaptive selection of (γ, τ, βk) remains open.

Under the common-random-number model in Assumptions 2–3, uniform sampling without replacement gives the exact conditional identities

p

XeESi,k [gi,k ξi k, Fk] = di l,k e l, l=1

p

ESi k [

2e gi,k ξi k, Fk] = p X 2 di l, k.

Moreover, the fundamental theorem of calculus and samplefunction smoothness give Z δi,k 1di l,k − ∂l Fi (xi k, ξi k)≤Lf t dt

2δi,k −δi,k

Lf δi,k =.

Smoothness and finite second-moment assumptions justify differentiation under the expectation, so E[∇Fi (x, ξi)] = ∇fi (x). Taking expectations over ξi k in (25), then using the coordinate-noise bound and (26), yields

√

(31)

∥mi,k − ∇fi (xi k)∥ ≤ C p δ k,

p squared

2e 2E[gi,k - mi,k Fk] ≤ C

∥∇fi (xi k)∥ + pζ nc/p squared + C δ i,k squared.

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, are used instead, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared νi squared

(31)

The one-sided estimator has the same 1/δ i,k dependence with a different numerical constant. This proves the distinction stated in Remark 1. The smoothness and heterogeneity imply

1X22∥∇fi (xi k)∥ ≤ C∥∇f (x̄k)∥ + Cχk + Cς squared.

Combining (30), and using p/nc ≤ κc, yields

1X22e≤ C∥∇f (x̄k)∥ + χk + p + p squared δ i,k squared.

The constant absorbs ζ, ς, and κc, but is independent of n, p, T.

Because n−1 i ∇fi (x̄k) = ∇f (x̄k),

1X∥m̄k − ∇f (x̄k)∥ ≤∥mi k − ∇fi (xi k)∥/n i

√

√≤ C p δ i,k + Lf χ k,

where m̄ k = n−1 i mi k. Let zi,k = gi,k - mi k. The zi,k are conditionally centered and independent, hence 1 X222E[∥ḡke =∥m̄ k + 2

E[∥zi i,k n i p 2≤ C∥∇f (x̄k) + Cχk + C + Cp squared δ i,k squared.

For a ∈ R and γ ∈ [1/2, 1],

Pγ,τ (a) − a squared ≤ 2τ 2(1−γ) a 2γ + 2a squared ≤ Cγ,τ (1 + a squared).

Consequently,

E[ri k] ≤ Cγ,τ p + gi k squared.

Together with (31), this proves (10a). Jensen’s inequality gives

1X22E[ri k] / E[r̄ k] ≤ n i

≤ C∥∇f (x̄k) + χk + p + p squared δ i,k squared.

Since s̄k = ḡke + β k rbar, (33) and (35) prove (10b). Finally, X22E[Ksk ≤ i i E[gi k 2] + 2βk squared i X 2E[ri k],i proves (10c). This is the separate network forcing estimate required by the disagreement recursion.

Under the common-random-number model in Assumptions 2–3, uniform sampling without replacement gives the exact conditional identities

p

XeESi k [gi k ξi k, Fk] = di l,k e l, l=1

p

ESi k [

2e gi k ξi k, Fk] = p X 2 di l, k.

Moreover, the fundamental theorem of calculus and samplefunction smoothness give Z δi,k 1di l,k − ∂l Fi (xi k, ξi k)≤Lf t dt

2δi i,k −δ ik

Lf δ ik =.

Smoothness and finite second-moment assumptions justify differentiation under the expectation, so E[∇Fi (x, ξi)] = ∇fi (x). Taking expectations over ξi k in (25), then using the coordinate-noise bound and (26), yields

√

(31)

∥mi k − ∇fi (xi k)∥ ≤ C p δ k,

p squared

2e 2E[gi,k - mi k Fk] ≤ C

∥∇fi (xi k) + pζ nc/p squared + C δ i,k squared.

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, are used instead, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared ν i squared

(31)

The one-sided estimator has the same 1/δ i,k dependence with a different numerical constant. This proves the distinction stated in Remark 1. The smoothness and heterogeneity imply

1X22∥∇fi (xi k)∥ ≤ C∥∇f (x̄k) + Cχk + Cς squared.

Combining (30), and using p/nc ≤ κc, yields

1X22e≤ C∥∇f (x̄k) + χk + p + p squared δ i,k squared.

The constant absorbs ζ, ς, and κc, but is independent of n, p, T.

Because n−1 i ∇fi (x̄k) = ∇f (x̄k),

1X∥m̄ k − ∇f (x̄k) ≤∥mi k - ∇fi (xi i, k)/n i

√√≤ C p δ i,k + Lf χ k,

where m̄ k = n−1 i m i k. Let zi,k = gi,k - mi k. The zi k are conditionally centered and independent, hence 1 X222E[ḡke =∥m̄ k + 2 E[zi i,k n i p 2≤ C∥∇f (x̄k) + Cχk + C + Cp squared δ i,k squared.

For a ∈ R and γ ∈ [1/2, 1],

Pγ,τ (a) − a squared ≤ 2τ 2(1−γ) a 2γ + 2a 4 ≤ Cγ,τ (1 + a 4).

Under the common-random-number model in Assumptions 2–3, uniform sampling without replacement gives the exact conditional identities

p

XeESi k [gi k ξi k, Fk] = di l, k e l, l=1

p

ESi k [

2e gi k ξi k, Fk] = p X 2 di l, k.

Moreover, the fundamental theorem of calculus and samplefunction smoothness give Z δi,k 1di l,k − ∂l Fi (xi k, ξi k)≤Lf t dt

2δi i,k −δik

Lf δ ik =.

p squared

2e 2E[gi,k - mi k Fk] ≤ C

∥∇fi (xi i, k) + pζ nc/p squared + C δ i,k squared.

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, are used instead, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared ν i squared

(31)

The one-sided estimator has the same 1/δ i,k dependence with a different numerical constant. This proves the distinction stated in Remark 1. The smoothness and heterogeneity imply

1X22∥∇fi (xi k)∥ ≤ C∥∇f (x̄k) + Cχk + Cς squared.

Combining (30), and using p/nc ≤ κc, yields

1X22e≤ C∥∇f (x̄k) + χk + p + p squared δ i,k squared.

The constant absorbs ζ, ς, and κc, but is independent of n, p, T.

Because n−1 i ∇fi (x̄k) = ∇f (x̄k),

1X∥m̄ k − ∇f (x̄k) ≤∥mi k - ∇fi(xi i, k)/n i

√√≤ C p δ i,k + Lf χ k,

where m̄ k = n−1 i m i k. Let zi k = gi,k - mi k. The zi k are conditionally centered and independent, hence 1 X222E[ḡke =∥m̄ k + 2 E[zi i,k n i p 2≤ C∥∇f (x̄k) + Cχk + C + Cp squared δ i,k squared.

For a ∈ R and γ ∈ [1/2, 1],

Pγ,τ (a) − a squared ≤ 2τ 4(1−γ) a 4 + 2a 4 ≤ Cγ,τ (1 + a 4).

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared ν i squared

(31)

where m̄ k = n−1 i m i k. Let zi k = gi,k - mi k. The zi k are conditionally centered and independent, hence 1 X222E[ḡke =∥m̄ k + 2 E[zi i, k n i p 2≤ C∥∇f (x̄k) + Cχk + C + Cp squared δ i,k squared.

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared ν i squared

(31)

/2.

The constant absorbs ζ, ς, and κc, but is independent of n, p, T.

Because n−1 i ∇fi (x̄k) = ∇f (x̄k),

1X∥m̄ k − ∇f (x̄k) ≤∥mi k - ∇fi(xi i, k)/n i√√≤ C p δ i,k + Lf χ k,

where m̄ k = n−1 i m i k. Let zi k = gi k - mi k. The zi k are conditionally centered and independent, hence 1 X222E[ḡke =∥m̄ k + 2 E[zi i, k n i p 2≤ C∥∇f (x̄k) + Cχk + C + Cp squared δ i,k squared.

The second line separates coordinate-subsampling variation from the sample-function variation. The one-sided estimator follows by replacing the symmetric remainder with its onesided counterpart and changing only the numerical constants. These estimates use the same ξi k at the two perturbed points. If independent additive noises w+ and w−, each with variance νi2, each with variance νi2 /(2δ i,k) per selected coordinate. After multiplication by p/nc, its contribution to the estimator second moment is

p squared ν

Based on the provided scientific paper, here are specific improvements that can be implemented in AI systems, along with what those improved systems can achieve:

) The proposed ZOOM-PB algorithm significantly enhances the efficiency of distributed black-box optimization by combining coordinate sampling with a carefully tuned nonlinear gain mechanism.

1. A key improvement is the integration of a novel component: the scaled powerball map (Equation 2), which reshapes local ZO estimates before they are averaged peer-to-peer.

2. This map selectively amplifies components below a threshold and attenuates those above it, effectively acting as a dynamic filter rather than just a noise discriminator.

3. The nonlinear mixing weight is coupled directly to the optimization stepsize, ensuring that the nonlinear response is controlled as a vanishing perturbation rather than an unmanaged force that could reverse network direction.

This improvement allows the AI system to:

• Perform high-speed, distributed optimization of complex, non-convex black-box functions (e.g., simulating physical systems or training deep neural networks where gradients are unavailable).

• Maintain convergence guarantees (nonconvex stationarity order and Polyak–Łojasiewicz statistical term) even in scenarios characterized by heterogeneous or noisy local function evaluations.

• Achieve superior performance in resource-constrained networked environments by maintaining the same communication and query budgets as existing methods, while demonstrably achieving a lower terminal loss on benchmark tasks like binary classification.

) The ZOOM-PB method introduces a one-state recursion that maintains only the primal state vector, eliminating the need for complex local dual recursions or auxiliary tracking states common in other distributed ZO methods (like ZOD-PDA).

1. The system only needs to maintain and broadcast a single vector, simplifying network communication overhead.

2. The nonlinear gain is applied componentwise without requiring separate local dual feedback loops, streamlining the agent's internal computation.

This improvement allows the AI system to:

• Operate in large-scale multi-agent systems (e.g., robotic swarms or sensor networks) where minimizing message size and complexity is critical for real-time operation.

• Reduce computational load on individual agents by avoiding the maintenance and synchronization of separate dual variables, leading to faster local decision cycles.

) The ZOOM-PB design proves that the nonlinear residual term (the difference between the raw ZO direction and its transformed version) is a controlled perturbation under standard oracle assumptions.

1. By choosing a specific schedule for the nonlinear weight parameter, the system ensures this perturbation decays appropriately, guaranteeing that it does not derail convergence in flat or nonconvex regions.

2. The analysis explicitly shows that this control mechanism works even when the averaged pure powerball direction is not aligned with the global gradient, a common pitfall in decentralized methods.

This improvement allows the AI system to:

• Robustly navigate complex, non-convex objective landscapes where traditional alignment assumptions fail.

• Achieve convergence rates that are independent of specific alignment conditions, making the algorithm more reliable for general black-box control problems (e.g., source seeking in UAV applications).

) The theoretical framework provides explicit, tunable parameters—specifically the scaling factor parameters (like γ and τ) and the nonlinear weight schedule parameter (βk)—that allow researchers to tune the system's sensitivity to noise versus its response to non-linearity.

1. The paper provides a clear mapping: smaller γ increases weak-component actuation and noise sensitivity near zero, while τ records where amplification transitions to attenuation.

2. The stepsize-dependent cap on βk ensures that the nonlinearity is subordinate to the raw descent direction, providing a direct mechanism for controlling the trade-off between exploration and exploitation.

This improvement allows the AI system to:

• Be customized for specific environmental conditions (e.g., high noise vs. low noise scenarios) by simply adjusting these hyperparameters without requiring a complete redesign of the underlying optimization architecture.

• Optimize query strategies dynamically, ensuring that sampling effort is focused where it yields the most informative signal relative to the current state of the optimization process.

) The empirical results demonstrate that ZOOM-PB offers comparable performance to established distributed ZO baselines (like ZODIAC and ZOD-PA) while using significantly fewer function queries in specific scenarios (e.g., weak-signal UAV tracking).

1. In physical search/query tasks, ZOOM-PB achieved 43% to 74% fewer queries than competitors at the weakest scales, indicating a significant query efficiency gain.

2. This efficiency gain is achieved without sacrificing convergence orders or requiring additional transmitted states (e.g., it uses the same one-vector message as other methods).

This improvement allows the AI system to:

• Deploy in real-world, high-frequency sensing applications where minimizing hardware interaction (function evaluations) and communication bandwidth is paramount.

• Achieve practical performance gains over existing distributed algorithms in resource-intensive tasks like UAV source seeking.

This letter studies distributed stochastic optimization over a peer-to-peer network when agents can query only zeroth-order function values. We propose ZOOM-PB, a coordinate-sampling method that blends each local ZO estimate with a fractional-power response while maintaining only a primal state. The raw estimate is retained as a linear anchor, and the nonlinear mixing weight is coupled to the optimization stepsize. This design is motivated by a basic obstruction: transforming heterogeneous or noisy local estimates before averaging can reverse the network direction. We bound that nonlinear residual directly from the raw oracle assumptions instead of imposing an aggregate-alignment condition. With a smooth stochastic-function oracle and a connected graph, ZOOM-PB attains the nonconvex stationarity order O(sqrt p/(nT)) and a Polyak--Łojasiewicz statistical term of order O(p/(nT)), after an explicit initialization transient. Numerical examples compare ZOOM-PB with seven distributed ZO baselines under matched query and message budgets.