Stochastic Mean-Field LQ Stackelberg Differential Games with Random Coefficients: Theory and a Deep FBSDE Picard Solver

Abstract: This paper studies a stochastic mean-field linear-quadratic Stackelberg differential game with random coefficients. The interaction between mean-field terms and random coefficients precludes the direct use of conventional decoupling techniques. We apply an extended Lagrange multiplier method to derive an affine operator representation of the follower's optimal response. The induced leader problem is then formulated as a generalized stochastic LQ control problem with operator-valued coefficients, and the Stackelberg optimal control is characterized through a Riccati-free coupled FBSDE system. We further develop a Deep FBSDE Picard Solver that preserves the Stackelberg order through follower-response learning, response-sensitivity extraction, leader optimization, and neural augmented Lagrangian enforcement of mean-field consistency constraints. Numerical studies covering convergence diagnostics, discretization sensitivity, Riccati calibration, ablation tests, local equilibrium validation, Stackelberg--Nash comparisons, and a financial application support the effectiveness of the proposed framework.

• The paper develops a Riccati-free extended Lagrange multiplier framework that represents the follower’s optimal control as an affine operator of the initial state and leader’s control, enabling a tractable leader problem despite random coefficients.
• The Deep FBSDE Picard Solver uses phase-separated neural networks, automatic differentiation, and augmented-Lagrangian updates to enforce mean-field consistency, achieving residuals as low as 2×10⁻⁴ and violations below 0.02.
• Numerical tests show that preserving Stackelberg sensitivities is essential, as removing them increases the leader’s cost by 49.3%, while unilateral-deviation tests support equilibrium accuracy within approximately 1%.

This paper studies a stochastic mean-field linear-quadratic (LQ) Stackelberg differential game with random coefficients and develops a deep-learning-based numerical solver for the resulting optimality system (2605.12950). The work combines two previously separate lines of research: the leader–follower stochastic LQ framework of Yong, in which the follower's rational response turns the leader's problem into an FBSDE-constrained control problem, and the extended Lagrange multiplier (ELM) method of Xiong and Xu for mean-field LQ control with random coefficients. The central technical obstacle is that random coefficients destroy the classical decoupling machinery: adjoint equations contain cross-moment terms such as $\mathbb{E}[A(t)^\top Y(t)]$ that cannot be factorized into products of expectations.

Problem formulation

On a filtered probability space with a one-dimensional Brownian motion, the state obeys a linear SDE whose drift and diffusion depend on $X(t)$, $\mathbb{E}[X(t)]$, both controls $u_1$ (follower) and $u_2$ (leader), and adapted random coefficient matrices $A_i, B_i, C_i, D_i$. Each player minimizes a quadratic cost with random weights $Q_i, R_i, G_i$ and their mean-field counterparts $\bar Q_i, \bar R_i$, over square-integrable adapted controls. Under assumptions (H1)–(H2)—boundedness of the operator-valued coefficients and uniform positivity of $R_i, \bar R_i$ together with $\bar Q_i \ge \delta I$—the paper establishes standard a priori estimates for the state equation and well-posedness of the associated linear BSDEs.

The Stackelberg structure is posed as two nested problems: (MFSOLQ-F), where the follower minimizes $X(t)$0 given a fixed leader control, and (MFSOLQ-L), where the leader minimizes $X(t)$1 anticipating the follower's optimal response map.

Follower's problem and affine response representation

Strict convexity of $X(t)$2 in $X(t)$3 is established via an operator-theoretic argument: expressing the state as an affine function of $X(t)$4 through bounded linear operators, the quadratic component of the cost is bounded below by $X(t)$5, which yields existence and uniqueness of the follower's optimum. The stochastic maximum principle characterizes the optimum through a coupled FBSDE with stationarity condition

$X(t)$6

Because the coefficients are not independent of $X(t)$7, Riccati decoupling is unavailable. The ELM method instead introduces deterministic auxiliary processes $X(t)$8, $X(t)$9 as constraints, relaxes them with extended Lagrange multipliers $\mathbb{E}[X(t)]$0, and exploits a min–max duality equality so that optimization order over controls and multipliers is immaterial. Solving the resulting sequence of subproblems (F-1)–(F-3) produces linear operators $\mathbb{E}[X(t)]$1 such that all relevant processes are affine in $\mathbb{E}[X(t)]$2; the optimal multipliers solve a linear operator equation, and the final pair $\mathbb{E}[X(t)]$3 solves a coercive block operator equation guaranteed solvable by (H2). When the multiplier system matrix $\mathbb{E}[X(t)]$4 is non-invertible, the paper handles the degenerate case by showing that any kernel perturbation of the multipliers maps to zero under the control-space projection, so the affine representation survives regardless.

The main structural result is that the follower's optimal response admits the form

$\mathbb{E}[X(t)]$5

with bounded linear operators $\mathbb{E}[X(t)]$6 and an inhomogeneous term. This affine response-operator structure is what makes the subsequent algorithmic design possible: it gives the leader explicit access to how its control propagates through the follower's rational behavior.

Leader's problem

Substituting the follower response into the dynamics yields a generalized MFSLQ problem with aggregated operator-valued coefficients $\mathbb{E}[X(t)]$7, $\mathbb{E}[X(t)]$8, which remain bounded under (H1). Applying the same ELM machinery characterizes the leader's optimal control by a coupled FBSDE with feedback

$\mathbb{E}[X(t)]$9

together with an operator equation for the leader's optimal expectation pair $u_1$0. The full Stackelberg equilibrium is thus characterized Riccati-free, through a deeply coupled forward–backward system with endogenous operator-valued coefficients.

The Deep FBSDE Picard Solver

On the computational side, the paper proposes the Deep FBSDE Picard Solver (DFPS), which deliberately preserves the sequential Stackelberg order rather than treating the bilevel system simultaneously. Three network families are used: AdjointNets approximate $u_1$1 conditioned on time, state, and a context vector $u_1$2 encoding all model coefficients; MacroNets approximate the mean-field quantities $u_1$3 from $u_1$4 only; LambdaNets parameterize the Lagrange multipliers. Conditioning on $u_1$5 allows a single trained model to cover many coefficient realizations without retraining.

A key design point concerns mean-field consistency. A naive Monte Carlo plug-in treats $u_1$6 and $u_1$7 as exogenous batch statistics, which fails to enforce the fixed-point consistency between macroscopic variables and policy-induced trajectories. DFPS instead imposes consistency constraints $u_1$8, $u_1$9 through an augmented Lagrangian scheme with proximal dual updates on LambdaNets and adaptive penalty growth (factor 1.1 upon violation stagnation). Proposition 4.1 gives an asymptotic feasibility guarantee: under uniformly bounded dual-subproblem and LambdaNet approximation errors, the constraint residual satisfies $u_2$0, so violations vanish if penalties diverge. This guarantee is conditional on Assumption 4.1 (bounded inexact-update errors); no rigorous a posteriori error analysis covering neural approximation error or Picard contraction under random operator coefficients is provided, which the authors state explicitly.

The three-stage pipeline trains the follower under exploratory leader controls, freezes it, extracts the response sensitivities $u_2$1 and $u_2$2 by automatic differentiation, and then trains the leader against the aggregated follower-induced dynamics. An optional fully coupled refinement loop exists but was never activated because tolerances were met without it.

Numerical results

The experimental program covers convergence, discretization sensitivity, a Riccati calibration, ablations, equilibrium validation, and a financial application. Several quantitative findings stand out:

• Residuals and feasibility: follower and leader BSDE residuals reach $u_2$3 and $u_2$4 within five Picard iterations; the follower adjoint terminal mismatch is $u_2$5 relative; all four mean-field violations terminate below the tolerance $u_2$6 (final values between $u_2$7 and $u_2$8). Across three seeds, costs are stable ($u_2$9, $A_i, B_i, C_i, D_i$0).
• Temporal convergence: under constant coefficients, the follower cost approaches the Riccati reference $A_i, B_i, C_i, D_i$1 with $A_i, B_i, C_i, D_i$2 relative error at $A_i, B_i, C_i, D_i$3, and log-log self-convergence slopes of $A_i, B_i, C_i, D_i$4 (self) and $A_i, B_i, C_i, D_i$5 (vs. Riccati) are consistent with first-order Euler–Maruyama accuracy. Under random coefficients, self-convergence improves from $A_i, B_i, C_i, D_i$6 at $A_i, B_i, C_i, D_i$7 to $A_i, B_i, C_i, D_i$8 at $A_i, B_i, C_i, D_i$9.
• Scaling: trainable parameters grow as $Q_i, R_i, G_i$0 ($Q_i, R_i, G_i$1), contrasted with exponential finite-difference grid growth; wall-clock time is dominated by fixed GPU overhead in the tested regime. The authors caution this is profiling, not a large-$Q_i, R_i, G_i$2 convergence analysis.
• Riccati sanity check: mean relative error of $Q_i, R_i, G_i$3 across three constant-coefficient scenarios, attributed to distributional shift since DFPS is trained on broad random ranges without targeted fine-tuning.
• Ablations: removing bilevel sensitivity ($Q_i, R_i, G_i$4, i.e., a Nash-type approximation) leaves $Q_i, R_i, G_i$5 unchanged but increases the leader's cost by $Q_i, R_i, G_i$6; simultaneous training without phase separation raises both costs ($Q_i, R_i, G_i$7, $Q_i, R_i, G_i$8) despite converging to a small BSDE loss of $Q_i, R_i, G_i$9—a notable demonstration that small residuals do not imply correct sequential structure; removing the augmented Lagrangian causes complete divergence.
• Equilibrium validation: unilateral-deviation tests with re-computed follower responses show cost increments within a $\bar Q_i, \bar R_i$0 tolerance band, with maxima of $\bar Q_i, \bar R_i$1 ($\bar Q_i, \bar R_i$2) and $\bar Q_i, \bar R_i$3 ($\bar Q_i, \bar R_i$4), supporting local Stackelberg optimality up to neural approximation accuracy.
• Financial application: in a mean-variance portfolio game (manager as leader, investor as follower) with stochastic volatility, the investor's tracking cost rises about $\bar Q_i, \bar R_i$5 relative to the deterministic baseline, while the manager's cost stays within Monte Carlo noise; across volatility regimes the differences are within one standard deviation and are presented as qualitative trends only. Notably, the presence of volatility uncertainty rather than its magnitude drives the investor's cost increase.

Limitations and open questions

Several limitations are conceded directly. The asymptotic feasibility result rests on Assumption 4.1, and the paper omits a posteriori error analysis encompassing neural approximation errors and Picard contraction under random operator-valued coefficients; a rigorous contraction justification would typically require a small-horizon condition that is not established here. The Riccati comparison carries a $\bar Q_i, \bar R_i$6 mean error attributable to distributional shift, and no closed-form baseline exists in the primary random-coefficient regime, so validation there relies on residual diagnostics and deviation tests rather than reference solutions. The financial experiment fixes coefficient instances per regime, and cross-regime cost differences are not statistically significant. The scalability study does not establish asymptotic convergence for high-dimensional states. Open questions include whether the sequential extraction procedure retains accuracy when the optional joint refinement becomes necessary, and how the framework behaves under genuinely time-varying adapted coefficients beyond the piecewise-constant realization used in experiments.

Conclusion

The paper provides a Riccati-free characterization of stochastic mean-field LQ Stackelberg equilibria with random coefficients via an extended Lagrange multiplier method, yielding affine operator representations of the follower's response and a coupled FBSDE system for the leader. The DFPS algorithm translates this structure into a phase-separated neural solver with augmented-Lagrangian enforcement of mean-field consistency, supported by an asymptotic feasibility proposition and extensive empirical validation including ablations and unilateral-deviation tests. The combination of theoretical response-operator structure and sensitivity extraction via automatic differentiation is the mechanism that lets the leader anticipate follower behavior without solving higher-order variational adjoints, and the reported results indicate the approach is effective in regimes where Riccati methods do not apply.