arXiv CSopen access

PAC-Bayesian Certificates for Quadratic Closed-Loop Control

arXiv:2606.28281v1 [eess.SY] 26 Jun 2026

PAC-Bayesian Certificates for Quadratic Closed-Loop Control Domagoj Herceg∗ Department of Mechanical Engineering Eindhoven University of Technology The Netherlands [email protected]

Abstract PAC-Bayesian bounds provide finite-sample guarantees for data-dependent randomized predictors, but applying them to learning-based control is difficult because the natural objective is a quadratic trajectory cost. Such losses are unbounded, non-Lipschitz , and lead to response-dependent Chernoff terms. We employ System Level Synthesis parameterization, which exposes the closed-loop trajectory map of a linear system directly and makes the quadratic control loss amenable to explicit certification. Moreover, we provide a set of PAC-Bayes-Chernoff certificates for posterior distributions over feasible closed-loop responses. For Gaussian disturbance trajectories with arbitrary covariance, we derive an exact one-sided Gaussian transform and a tractable quadratic upper bound expressed through closed-loop sensitivity quantities. We also derive a posterior-localized surrogate for settings where pointwise closed-loop response certificates are unavailable or have support related admissibility issues. Although PAC-Bayes certifies a non-degenerate posterior, the convex quadratic form of the SLS loss transfers the certificate to the posterior mean response. We present a deterministic mean response deployment result that is particularly suitable for control while retaining the stochastic posterior in the bound. Additionally, we provide a data-driven bound for this deployment, transitioning away from an oracle bound. Minimizing this bound naturally results in a learning algorithm for control selection from data. Numerical experiments on a double integrator show that the algorithm acts as a sensitivity-aware finite-sample regularizer, improving held-out cost and reducing closed-loop sensitivity in the low-data regime.

1

Introduction

Learning methods have made great strides in the last decade, mostly due to the enormous availability of data. Due to the complexity of modern control problems, there has been notable interest in datadriven/learning methods in the control community [Verheijen et al., 2023, Dörfler et al., 2023, Recht, 2019, Markovsky and Dörfler, 2021]. Unlike the usual bounded costs, such as zero-one training loss in machine learning, squared unbounded costs are the norm in control applications for the most part. This stems from a connection to stability via Lyapunov quadratic forms [Grüne and Pannek, 2016], the explicit solution for the canonical Linear Quadratic Regulator problem [Athans and Falb, 2013], and the general convexity of the cost [Boyd and Vandenberghe, 2004], which enables fast and reliable solution algorithms. We are often given an approximate model of the system, and the trajectory of the control actions is optimized w.r.t the loss induced by that model. In what follows, we will interlace the terminology from both control and ML, which may result in slightly unfamiliar notation for both ∗ Fallback address: [email protected]

Preprint.

sides. However, we lean more toward the machine learning notation borrowed from PAC-Bayesian literature. Classical linear control typically certifies a single controller, while robust control certifies against a worst-case set. In learning-based control, however, controllers are often selected from finite rollout data, and the pressing question is whether the selected controller generalizes well to unseen disturbance trajectories. As a consequence, it is important to offer guarantees that the learned controller will perform well on unseen data. In the Probably Approximately Correct (PAC) framework [Valiant, 1984], the goal is to provide that guarantee with high probability. In our view, a particular instance of the PAC framework, the PAC-Bayes theory[Alquier, 2021, Guedj, 2019], provides a natural intermediate guarantee as it certifies the expected deployment cost of a posterior distribution over hypotheses, with a KL penalty measuring how far this posterior moves from a data-independent prior. This gives a finite-sample certificate for a distribution of hypotheses rather than a single point estimate and can be interpreted as a robustness measure regarding perturbations. Other PAC (or PAC style) methods usually found in control include conformal predictions [Lindemann et al., 2024, Angelopoulos et al., 2024, Vlahakis et al., 2024] and the scenario approach [Campi and Garatti, 2018, Campi et al., 2009]. Unlike the traditional view of designing a controller directly, System Level Synthesis [Wang et al., 2019] (SLS) makes this PAC-Bayesian viewpoint tractable for linear systems. By parameterizing feasible closed-loop responses directly, the weighted trajectory can be written as an affine function of the disturbance trajectory, so the finite-horizon quadratic control cost is exactly the squared norm of this response. This exposes the closed-loop sensitivity maps entering the certificate and allows explicit one-sided Chernoff calculations for Gaussian disturbance models. The SLS parameterization makes this certificate remarkably interpretable: posterior uncertainty is penalized in directions that strongly amplify disturbances into quadratic state-input costs and can remain spread in directions that are benign for closed-loop performance. Moreover, it certifies that the learned posterior over closed-loop responses generalizes from sampled disturbance trajectories to new trajectories drawn from the same deployment law. Related work Despite its remarkable success in the field of machine learning, the adoption of PACBayesian methods in learning control has remained somewhat muted. PAC-Bayesian methods have only recently been used to certify learning-based controllers by translating finite-sample generalization guarantees from supervised learning to control. Work by [Majumdar and Goldstein, 2018, Majumdar et al., 2020] introduced PAC-Bayes control for robotic policies that generalize across environments, with algorithms that optimize PAC-Bayes bounds over finite or continuously parameterized policy classes. More recent work [Boroujeni et al., 2025, 2024] has developed PAC-Bayesian optimal control frameworks for stochastic nonlinear systems, including posterior-based controller design and stabilizing neural controller parameterizations. These works show that PAC-Bayes is a natural tool for certifying data-dependent randomized controllers, but they primarily focus on bounded costs. The latter works also focus on stability-preserving nonlinear control, assuming a known stabilizing controller, in addition to considering bounded disturbances. Rather than targeting general nonlinear policy classes, we exploit the finite-horizon affine SLS parameterization for linear systems, which exposes the closed-loop trajectory as an affine function of the disturbance. This allows us to treat the native quadratic state-input cost directly. In particular, we derive exact and tractable PAC-Bayes-Chernoff certificates for unbounded quadratic losses under Gaussian disturbance trajectories, with arbitrary covariance.A related work to ours is a recent preprint [Herceg and Antunes, 2026] that combines PAC-Bayes with distributional robustness to address the mismatch between training and deployment environments via SLS parameterization. However, they focus on norm type losses, heavily exploit the Lipschitz nature of the loss, and use SLS as a means to an end without explaining the structural benefits of such a choice. In contrast, our focus is on the characterization of the true quadratic closed-loop loss under a fixed deployment disturbance law, which was highlighted as a challenge in their work. The Gaussian certificate developed here uses the one-sided quadratic transform rather than a Lipschitz or sub-Gaussian proxy, yielding sensitivity terms based on the induced closed-loop covariance. Additionally, we provide a data-driven algorithm transitioning away from oracle bounds and fully justify the use of the SLS framework. Contributions. We make four contributions: (i) We show that finite-horizon SLS provides a natural setting for the stochastic PAC-Bayes certification of quadratic closed-loop control losses. (ii) Building on the above, we derive exact and tractable one-sided PAC-Bayes-Chernoff certificates 2

for unbounded quadratic closed-loop losses under Gaussian disturbance trajectories with arbitrary covariance. (iii) We also prove a deterministic mean-response certificate. The non-degenerate posterior makes the complexity term well-defined, while the quadratic SLS identity transfers the randomized certificate to the deployed posterior mean, up to a curvature mismatch term (empirical to population mismatch), blending deterministic deployment with stochastic certification. (iv) We derive a data-driven certificate by inflating the empirical covariance. Hence, we turn an oracle bound into a computable one. Additionally, we show empirically on a double integrator that posterior optimization followed by mean deployment outperforms point optimization in the low to mid data regime, yielding lower held-out cost and lower closed-loop sensitivity.

2

PAC-Bayes Setup and Quadratic SLS Losses

In this section we introduce the necessary preliminaries regarding PAC-Bayes theory and the SLS parameterization. Moreover, we discuss how the SLS parametrization naturally blends with PACBayes, a connection not made before in the literature. 2.1

PAC-Bayes

PAC-Bayesian methods [Catoni, 2007, McAllester, 1998] have garnered a lot of attention due to their powerful ability to provide generalization guarantees for randomized predictors, with the most notable showcase being the first non-vacuous bounds for deep neural networks [Dziugaite and Roy, 2017, Pérez-Ortiz et al., 2021]. A typical PAC-Bayesian inequality aims to bound the population risk, the performance after deployment, by a sum of expected empirical performance and a complexity term. Let us denote the space of all hypothesis/control policies/closed-loop responses with Θ ⊆ Rd . We will emphasize this overloaded nomenclature often throughout the paper to highlight the bridge between the disciplines. Our goal is to infer a distribution ρ ∈ M1 (Θ), called the posterior distribution. We briefly note that the posterior need not be the Bayesian posterior, see [Germain et al., 2016]. This inference is based on a dataset S = {wi }ni=1 where the data is assumed to be i.i.d. sampled from an unknown data-generating distribution D supported on W ⊆ Rw . Then the abstract loss function can be stated as ℓ : Θ × W 7→ R+ .

(1)

We define the population risk and the empirical risk (computable from data) as n 1X b ℓ(θ, wi ). L(θ) := n i=1

L(θ) := Ew∼D [ℓ(θ, w)]

(2)

PAC-Bayes inequality then bounds the difference between the two after taking the expectation over b the posterior θ ∼ ρ. To ease the notation, we will often write Eρ [L(θ)] and Eρ [L(θ)] with the implicit understanding that we are averaging over the variable in Θ, or θ ∼ ρ in the above. In addition, we need to specify a data-independent prior distribution π ∈ M1 (Θ) that will act as an anchor against overfitting by penalizing the deviation of the posterior w.r.t the prior. Finally, we specify the confidence δ ∈ (0, 1], indicating that the bound holds with a probability of no less than 1 − δ. PAC-Bayes Chernoff bounds There are different types of PAC-Bayes bounds, Langford-SeegerMaurer type [Langford and Seeger, 2001, Maurer, 2004], McAllester style McAllester [2003], and the literature on the subject is vast. We focus on typically more modern results for unbounded costs and, more specifically, on the results in [Casado et al., 2024]. The core idea of their approach is to bound the cumulant generating function (CGF) by a function ψ. Assumption 1 (θ-dependent bounded CGF). For every θ ∈ Θ, define the centered CGF 

 Λθ (λ) := log Ew∼D exp λ L(θ) − ℓ(θ, w) .

(3)

There exist b > 0 and a nonnegative convex function ψ(θ, λ) such that ψ(θ, 0) = ψ ′ (θ, 0) = 0 and Λθ (λ) ≤ ψ(θ, λ),

∀λ ∈ [0, b).

Assumption 1 is essential for the derivation of the following theorem. 3

(4)

Theorem 2 (PAC-Bayes bound under bounded CGF [Casado et al., 2024]). Suppose Assumption 1 holds. Then for any δ ∈ (0, 1), with probability at least 1 − δ over the draw of S ∼ Dn , the following holds simultaneously for all posteriors ρ ≪ π:   KL(ρ∥π) + log(n/δ) Eθ∼ρ [ψ(θ, λ)] b Eθ∼ρ [L(θ)] ≤ Eθ∼ρ [L(θ)] + inf + , (5) λ(n − 1) λ λ∈[0,b) | {z } Comp(ρ,S,δ)

There are a few things to note here that distinguish (5) from a more typical PAC Bayes bound of this type. First, the minimization in λ is exact and avoids any clumsy union bound arguments over a predetermined grid. Even more crucially, this bound avoids taking the expectation of the exponential of ψ under the prior and directly uses the posterior ρ in the term Eθ∼ρ [ψ(θ, λ)]. This is particularly well suited for control, as the cost will depend on the hypothesis/controller/closed-map distribution and not on the prior. Hence, we obtain a bound that depends on the performance of the actual optimized posterior we deploy. The term under the inf operator is the complexity term in which we can recognize KL(ρ∥π), an information-centric measure that punishes the information gain of the posterior with respect to the prior. 2.2

System Level Synthesis

To enable a tractable closed-loop and learning-based analysis, we adopt a finite-horizon System Level Synthesis (SLS) [Wang et al., 2019] parameterization for control. The key idea is to directly parameterize the closed-loop response from disturbances to state and input trajectories instead of the standard mapping through control policies. The original SLS formulation was developed for linear state feedback policies, which were later extended to affine policies [Sieber et al., 2021, Schüepp et al., 2025]. The motivation was the well known fact that quadratic Model Predictive Control (MPC) [Rawlings et al., 2020] with polyhedral constraints admits an affine feedback law depending on the region of the state space [Bemporad et al., 2002]. For more details about SLS parametrization, see Appendix B. Linear dynamical system

Consider a linear time invariant (LTI) system over a finite horizon

x̄t+1 = Ax̄t + B ūt + w̄t ,

t = 0, . . . , T − 1.

(6)

with x̄t ∈ Rnx , ūt ∈ Rnu , w̄t ∈ Rnx . The standard control approach would be to parameterize the controller as state (static for simplicity) feedback ūt = K x̄t , K ∈ Rnu ×nx , and rewrite the dynamics as x̄t+1 = (A + BK)x̄t + w̄t . However, this makes the bound less explicit. Even though it seems natural, this introduces a non-convex landscape concerning system design, as trajectories depend on polynomials of (A + BK). SLS parametrization for the stacked dynamics The evolution of system (6) can be written more compactly using stacked notation. Define the stacked vectors ⊤

x := [x̄0 x1 . . . x̄T ] ,

⊤

u := [ū0 u1 · · · ūT −1 ] ,

⊤

w := [x̄0 w̄0 · · · w̄T −1 ] .

of appropriate dimensions. Each sample wi is a finite-horizon disturbance trajectory. We note that the coordinates of a trajectory need not be independent. Then, the system closed-loop maps from disturbance w to the state and input trajectories can be written as          x Φx ϕ Φ x ϕx w = w+ x = (7) u Φu ϕu Φu ϕ u 1 where Φx ∈ R(T +1)nx ×(T +1)nx and Φu ∈ RT nu ×(T +1)nx are block lower-triangular operators encoding causality. Affine terms are ϕx ∈ R(T +1)nx , ϕu ∈ RT nu . Instead of solving for the control policies , we are interested in optimizing directly over the closed-loop response map represented by a tuple θ := {Φx , Φu , ϕx , ϕu }. However, response maps cannot be freely chosen as desired. Any closed-loop response tuple must obey the dynamics imposed by (6). More precisely, the closed-loop response tuple must be consistent with the linear achievability constraint   Φ ϕx [(I − ZA) −ZB] x = [I d] , (8) Φu ϕ u 4

where Z is the block-downshift operator, and A, B denotes the block-diagonal lifted system matrices. In our setting d = 0, but in general, it need not be so. What is important is that (8) is a linear constraint in the tuple, and it parameterizes all closed-loop responses. Moreover, the linear part of a time-varying control policy can be recovered from the closed loop responses, see Appendix B for further explanations and definitions of Z, A and B. 2.3

System Level Synthesis for Stochastic Quadratic Control

In this subsection, we take a step further than the standard literature and describe how SLS is the natural coordinate system for stochastic deployment in the native quadratic cost setting. To the best of our knowledge, we are the first to do so. Nullspace coordinates and feasible response directions. Vectorize the response tuple as   vec(Φx ) vec(Φu ) ϑ :=  . φx  φu

(9)

After incorporating the causal sparsity constraints, the achievability condition (8) can be written as a finite dimensional linear system Hϑ = h. This was already noted in [Herceg and Antunes, 2026] and the authors stopped there as their goal was to have a basis that can be used with posteriors that have unbounded support. Here, instead, we dive deeper into the particular properties of the quadratic cost and offer a much more complete characterization. Let ϑ0 be one feasible solution of Hϑ = h, and let N = [N1

···

N2

Np ]

(10)

be a fixed basis for ker(H). Then every feasible causal affine response can be written as ϑ(α) = ϑ0 + N α = ϑ0 +

p X

α ∈ Rp .

αk Nk ,

(11)

k=1

Pp Indeed, Hϑ(α) = Hϑ0 + k=1 αk HNk = h. Thus, every α ∈ Rp corresponds to a feasible causal closed-loop response. We place the PAC-Bayes prior and posterior over these free coordinates. The basis N fixes the coordinate system in which prior and posterior complexity are measured and is chosen before observing the data. Further, each nullspace direction Nk can be partitioned as   [k]  vec Φx   [k]  vec Φu  , Nk =  k = 1, . . . , p, (12)    φ[k]  x [k] φu [0]

[0]

[0]

[0]

while ϑ0 determines the corresponding objects Φx , Φu , φx , φu and Φx (α) = Φ[0] x + φx (α) = φ[0] x +

p X k=1 p X

αk Φ[k] x ,

Φu (α) = Φ[0] u +

αk φ[k] x ,

φu (α) = φ[0] u +

k=1

p X k=1 p X

αk Φ[k] u , (13) αk φ[k] u .

k=1

Weighted response maps and the two affine structures. Q̄ := blkdiag(Q0 , . . . , P ),

Let

R̄ := blkdiag(R0 , . . . , RT −1 ),

where P may be the terminal weight. Define the weighted closed-loop trajectory  1/2  Q̄ x yα (w) := = M (α)w + m(α), R̄1/2 u where M (α) :=

 1/2  Q̄ Φx (α) , R̄1/2 Φu (α)

m(α) :=

5

 1/2  Q̄ φx (α) . R̄1/2 φu (α)

(14)

(15)

The quadratic finite-horizon control loss is therefore ℓ(α, w) := ∥yα (w)∥22 = ∥M (α)w + m(α)∥22 .

(16)

By (13), the weighted response maps are affine in α: M (α) = M0 +

p X

αk Mk ,

m(α) = m0 +

k=1

where

"

# [k] Q̄1/2 Φx Mk := [k] , R̄1/2 Φu

p X

αk mk ,

(17)

k = 0, . . . , p.

(18)

k=1

"

# [k] Q̄1/2 φx mk := [k] , R̄1/2 φu

The pair (Mk , mk ) has a direct closed-loop interpretation as it is the cost weighted change in the disturbance induced map in the k-th feasible nullspace direction Nk . The representation above gives two complementary affine structures. First, for fixed α, the map w 7→ M (α)w + m(α) is affine in the disturbance trajectory. Therefore, ℓ(α, w) is quadratic in w. Under Gaussian disturbances, this is precisely the structure used in Section 4.1 to derive an exact one-sided Gaussian Chernoff transform. Second, for fixed w, (17) gives M (α)w + m(α) = M0 w + m0 +

p X

αk (Mk w + mk ) = a0 (w) + A(w)α,

(19)

···

(20)

k=1

where a0 (w) := M0 w + m0 ,

A(w) := [M1 w + m1

Mp w + mp ] .

Consequently, ℓ(α, w) = ∥a0 (w) + A(w)α∥22

(21)

is convex quadratic in α. Hence, both the population and empirical risks are convex quadratic functions of the SLS coordinate. This second structure is used in Section 5 to transfer a randomized PAC-Bayes certificate to the deterministic mean response. Takeaway In summary, SLS is central to our approach because it simultaneously guarantees multiple desirable properties: (i) the nullspace parameterization guarantees closed-loop feasibility for posterior distributions with unbounded support, (ii) it makes the loss quadratic in the disturbance trajectory for each fixed response, which enables the certificate in Section 4.1, (iii) it makes the loss quadratic in the free response coordinate for each fixed disturbance realization, which enables deterministic deployment in Section 5 , and finally (iv) the feasible basis directions (Mk , mk ) describe how posterior uncertainty over response coordinates maps the disturbances to cost, which is used in the data-driven certificate to quantify the additional covariance-mismatch penalty (see Section 5).

3

Posterior-localized PAC-Bayes bounds

Theorem 2 is directly useful when one has a pointwise certificate Λα (λ) ≤ ψ(α, λ) valid on a common interval [0, b). If the posterior distribution has unbounded support, the common admissible range over the posterior support may collapse even when most posterior mass is concentrated on "good" responses. For sub-Gaussian disturbances, Hanson-Wright [Hanson and Wright, 1971] yields such a certificate but only on a posterior-dependent interval λ ∈ [0, (Cσ 2 ∥M (α)∥2op )−1 ), which may collapse to zero under broad unbounded posteriors whenever the posterior assigns mass to directions along which ∥M (α)∥op is unbounded. Corollary 3 targets exactly this setting. 3.1

A posterior-localized Chernoff surrogate

Instead of first certifying each response and then averaging the certificate, we upper-bound the posterior average of the response-wise CGF by a single joint exponential moment under the deployed posterior. 6

Corollary 3 (Posterior-localized Chernoff surrogate). Let S = (w1 , . . . , wn ) ∼ Dn with n ≥ 2, and suppose that ℓ(α, w) ≥ 0 and L(α) < ∞ for every α ∈ Θ. Define n 1X b LS (α) := ℓ(α, wi ), n i=1

and, for λ ≥ 0, 

 Λα (λ) := log Ew∼D exp λ L(α) − ℓ(α, w) . For a posterior ρ ≪ π, define the joint posterior-localized CGF 

 b ρ (λ) := log Eα∼ρ, w∼D exp λ L(α) − ℓ(α, w) Λ , and its effective domain n o b ρ (λ) < ∞ . Ibρ := λ > 0 : Λ Then, with a probability of at least 1 − δ over S ∼ Dn , simultaneously for all ρ ≪ π with Ibρ ̸= ∅, ) ( b ρ (λ) Λ KL(ρ∥π) + log(n/δ) b S (α)] + inf + . Eα∼ρ [L(α)] ≤ Eα∼ρ [L (n − 1)λ λ λ∈Ibρ Proof. Since ℓ(α, w) ≥ 0, for every α ∈ Θ and λ ≥ 0, 
 Λα (λ) ≤ log Ew∼D exp λL(α) = λL(α) < ∞. Moreover, Λα is convex, nonnegative, and satisfies Λα (0) = Λ′α (0) = 0. Hence, Assumption 1 holds with b = ∞ by trivially choosing ψ(α, λ) = Λα (λ). Applying Theorem 2 yields, with probability at least 1 − δ, simultaneously for all ρ ≪ π,   KL(ρ∥π) + log(n/δ) Eα∼ρ [Λα (λ)] b Eα∼ρ [L(α)] ≤ Eα∼ρ [LS (α)] + inf + . λ>0 (n − 1)λ λ b ρ , and using Tonelli [Schilling, 2017] to exchange Now fix λ ∈ Ibρ . By the definitions of Λα and Λ expectations, we obtain 
 


b ρ (λ) < ∞. Eα∼ρ exp Λα (λ) = Eα∼ρ, w∼D exp λ L(α) − ℓ(α, w) = exp Λ Therefore, Jensen’s inequality gives 
 b ρ (λ). Eα∼ρ [Λα (λ)] ≤ log Eα∼ρ exp Λα (λ) = Λ Combining this inequality with the bound from Theorem 2 gives, for every λ ∈ Ibρ , b b S (α)] + KL(ρ∥π) + log(n/δ) + Λρ (λ) . Eα∼ρ [L(α)] ≤ Eα∼ρ [L (n − 1)λ λ Taking the infimum over λ ∈ Ibρ proves the result. We call the bound posterior-localized because admissibility is imposed on the joint law (α, w) ∼ ρ × D, rather than pointwise over every response in the posterior support. When explicit pointwise certificates are available, the averaged-CGF route in Theorem 2 is sharper. The role of Corollary 3 is to provide a fallback when response certificates are unavailable or have support admissibility issues. Remark 4 (Posterior-level admissibility). The admissibility requirement in Corollary 3 is posteriorlocalized. In particular, it does not require a common Chernoff domain over the entire feasible SLS class. Rather, the bound is meaningful for those posteriors ρ for which Ibρ is nonempty. The feasible SLS class may contain responses with poor behavior, while the learned posterior may concentrate on a particularly desirable neighborhood of closed-loop responses. The joint posterior-CGFs form is naturally conservative relative to the posterior average of pointwise CGFs since it uses Jensen’s inequality. 7

Quadratic SLS specialization. If L(α) = α0 + 2g ⊤ α + α⊤ Gα, G ⪰ 0, and ρ = N (µρ , Σρ ), 1/2 1/2 then Eρ eλL(α) < ∞ whenever I − 2λΣρ GΣρ ≻ 0. This condition depends on the deployed posterior covariance, not on the prior. However, it is conservative because it uses L − ℓ ≤ L and is mainly used as an example to highlight how to avoid radius collapse. Sharper certificates are given in the next section.

4

Gaussian certificates

Corollary 3 is distribution-free at the PAC-Bayes level, but it is defined for values of λ for which the posterior-localized exponential moment is finite. To obtain explicit, computable certificates, we now specialize for the Gaussian disturbance, which yields an exact closed-form expression for Λα (λ) by exploiting the quadratic-Gaussian structure of the loss. A similar approach can be taken for bounded disturbances, which yield a distribution-free certificate via Hoeffding’s lemma [Boucheron et al., 2013], which we do not cover here. In both cases, the resulting bound is obtained by substituting the corresponding certificate into Theorem 2. 4.1

Gaussian disturbance certificates

Next, we derive an exact expression for the inner exponential moment under a Gaussian disturbance model. Unlike standard Hanson-Wright bounds [Rudelson and Vershynin, 2013, Ziemann, 2025], which are typically stated for quadratic forms of vectors with independent sub-Gaussian coordinates, our Gaussian certificate allows a full trajectory-level covariance Σw . Thus, the disturbance sequence may be temporally correlated, and the certificate depends on the induced closed-loop covariance Σy (α) = M (α)Σw M (α)⊤ . Proposition 5 (Exact Gaussian integrand). Fix α ∈ Θ and let w ∼ D = N (µw , Σw ), so that yα := M (α)w + m(α) ∼ N (µy (α), Σy (α)) with µy (α) = M (α)µw + m(α) and Σy (α) = M (α)Σw M (α)⊤ , hence ℓ(α, w) = ∥yα ∥22 , L(α) = tr(Σy (α)) + ∥µy (α)∥22 . Define a shorthand Σ(α) = I + 2λΣy (α), then for every λ ≥ 0, h i
−1/2 −λ µy (α)⊤ Σ(α)−1 µy (α)
e Ew∼D eλ(L(α)−ℓ(α,w)) = eλL(α) det Σ(α) . (22) Equivalently, the centered CGF is Λα (λ) = λL(α) − 21 log det(I + 2λΣy (α)) − λ µy (α)⊤ (I + 2λΣy (α))−1 µy (α).

(23)

Proof. Since yα = M (α)w + m(α) ∼ N (µy (α), Σy (α)), we can write ℓ(α, w) = ∥yα ∥22 ,

L(α) = tr(Σy (α)) + ∥µy (α)∥22 .

Hence

2

Ew eλ(L(α)−ℓ(α,w)) = eλL(α) Eyα e−λ∥yα ∥2 . For a Gaussian vector y ∼ N (µ, Σ) with Σ ≻ 0 and any λ ≥ 0,   2 Ee−λ∥y∥2 = det(I + 2λΣ)−1/2 exp − λ µ⊤ (I + 2λΣ)−1 µ by [Mathai and Provost, 1992, Corollary 3.2a.2], but the same formula holds for Σ ⪰ 0 due to I + 2λΣ being well defined, which is a consequence of the one-sided CGF we consider (see Appendix A.1). Applying this identity with (µ, Σ) = (µy (α), Σy (α)) yields (22). Taking logarithms gives (23). Proposition 5 gives an exact expression for the one-sided CGF of the squared affine-SLS loss under Gaussian disturbances. Its importance is twofold. First, it preserves the damping induced by the negative loss term −ℓ(α, w). Second, the resulting expression is finite for every λ ≥ 0, since it depends on I + 2λΣy (α). Thus, in the Gaussian disturbance case, the one-sided Chernoff transform avoids a posterior-dependent finite radius. The certificate depends on the mean and covariance of the weighted closed-loop trajectory, µy (α) and Σy (α), providing a direct closed-loop sensitivity interpretation. 8

Corollary 6 (Quadratic upper bound for the Gaussian CGF). Under the assumptions of Proposition 5, for every α ∈ Θ and every λ ≥ 0,
Λα (λ) ≤ λ2 ∥Σy (α)∥2F + 2µy (α)⊤ Σy (α)µy (α) =: ψ(α, λ) (24) where

µy (α) := M (α)µw + m(α), Σy (α) := M (α)Σw M (α)⊤ . and therefore Theorem 2 yields the bound  h i KL(ρ∥π) + log(n/δ) b Eρ [L(α)] ≤ Eρ [L(α)] + inf + λEρ ϕ(α) , λ>0 λ(n − 1)

(25)

where ϕ(α) := ∥Σy (α)∥2F + 2µy (α)⊤ Σy (α)µy (α). Proof. We sketch the proof here and provide the full one in Appendix A.2. We start from eq. (23) and group the terms with µy (α) and Σy (α) separately. For the mean term, we use the identity I − (I + 2λΣy (α))−1 = 2λΣy (α)(I + 2λΣy (α))−1 to obtain 2λ2 µy (α)⊤ Σy (α)(I + 2λΣy (α))−1 µy (α) ≤ 2λ2 µy (α)⊤ Σy (α)µy (α). In the above, we used the fact that (I + 2λΣy (α))−1 ⪯ I. For the covariance term, let s1 , . . . , sr denote the eigenvalues of Σy (α). Then  r  X 1 1 λsi − log(1 + 2λsi ) . (26) λ tr(Σy (α)) − log det(I + 2λΣy (α)) = 2 2 i=1 A valid scalar inequality for x ≥ 0 is x − 12 log(1 + 2x) ≤ x2 (see Appendix A.2). Substituting x = λsi , we can bound (26) by λ2 ∥Σy (α)∥2F , which follows from the definition of the Frobenius norm. Combining the two estimates proves (34). Moreover, ψ(α, 0) = ψ ′ (α, 0) = 0 and ψ(α, λ) is clearly convex in λ satisfying Assumption 1 with b = ∞. Both terms of the certificate are closed-loop sensitivity quantities expressed through the induced covariance Σy (α). The covariance term ∥Σy (α)∥2F penalizes the overall disturbance amplification, while the mean term µy (α)⊤ Σy (α)µy (α) penalizes the mean offset in the same induced metric.

5

Deterministic mean-response deployment and data-driven bounds

PAC-Bayes certifies randomized posterior risks, whereas control applications typically deploy a deterministic controller. This stems from stringent testing and repeatability analysis. In this section, we will therefore use the posterior ρ as a certification object during learning but deploy only its mean response µρ . First, note that the controller that achieves it is Kmean = Φu (µρ )Φx (µρ )−1 . Furthermore, it needs to be said that we will keep the posterior distribution in the complexity term, blending stochastic certification and deterministic deployment, keeping the best of both worlds (so to speak) in a single formulation. This also avoids applying PAC-Bayes directly to the Dirac posterior δµρ , which would generally have infinite KL divergence with respect to a continuous prior. In the affine SLS parameterization, M (α) and m(α) are affine in α. Hence, for every fixed disturbance trajectory w, the loss ℓ(α, w) = ∥M (α)w + m(α)∥22 is convex quadratic in α. Assuming ED ∥w∥2 < ∞, the population risk is L(α) = Ew∼D [ℓ(α, w)] = α⊤ Gα α + 2gα⊤ α + α0 ,

Gα ⪰ 0.

(27)

Similarly, the empirical risk has a quadratic form b S (α) = α⊤ G b S α + 2b L gS⊤ α + α b0,S ,

b S ⪰ 0. G

(28)

We can now state the following proposition. Proposition 7 (Deterministic mean-response certificate). Let ρ be a posterior with a finite second moment, and denote the mean as µρ and the covariance as Σρ . Suppose that, with probability at least b S (α)] + Comp(ρ, S, δ) holds. 1 − δ, the randomized PAC-Bayes certificate Eα∼ρ [L(α)] ≤ Eα∼ρ [L Then, on the same event,   b S (µρ ) + Comp(ρ, S, δ) + tr (G b S − G)Σρ . L(µρ ) ≤ L 9

b S (α)] = L b S (µρ ) + Proof. Since L is quadratic, Eα∼ρ [L(α)] = L(µρ ) + tr(GΣρ ). Likewise, Eα∼ρ [L b S Σρ ). Substituting both identities into the randomized PAC-Bayes certificate and rearranging tr(G gives the claim. Note that this bound is still an oracle bound, but the unknown term can be bounded via standard concentration inequalities [Vershynin, 2018, Wainwright, 2019]. Proposition 7 shows that posterior learning can be seen as a distributional regularization of deterministic empirical risk minimization b S (α)] = L b S (µρ ) + (ERM). During training, the empirical posterior risk decomposes as Eα∼ρ [L b S Σρ ), so the posterior regularizes the deployed mean by encouraging a finite KL neighborhood tr(G around it to have low empirical cost. After transferring the certificate to the deterministic mean response, the full empirical spread penalty is not paid as a deployment cost, it is only the empirical to population curvature mismatch along Σρ that contributes. In particular, if the empirical and b S − G)Σρ = 0, then L(µρ ) ≤ population curvatures agree along the posterior covariance, i.e. tr(G b S (µρ ) + Comp(ρ, S, δ). L Thus, the deployed controller is deterministic, but its certificate is inherited from the non-degenerate posterior surrounding it. This is the mechanism by which PAC-Bayes provides a finite-sample certificate for mean deployment without incurring the infinite KL cost of a Dirac posterior. Proposition 8 (Data-driven deterministic certificate). Consider the zero mean Gaussian disturbance i.i.d. bw − setting wi ∼ N (0, Σw ) with m(α) = 0, and suppose that the covariance event EΣ = {∥Σ Σw ∥op ≤ ϵΣ } holds. Define the data-driven sensitivity coefficient, analogous to Corollary 6, as   2 bρ := Eα∼ρ M (α)(Σ b w + ϵΣ I)M (α)⊤ C . F

Let (MG )ij = tr (Mi⊤ Mj ) be the Gram matrix. Then, on the intersection of the PAC-Bayes event and the covariance event, the deterministic mean response µρ = Eρ [α] satisfies s b b S (µρ ) + 2 Cρ (KL(ρ∥π) + log(n/δPB )) + ϵΣ tr(Σρ MG ). L(µρ ) ≤ L (29) n−1 Finally, if P(EΣ ) ≥ 1 − δΣ , then the bound (29) holds with a probability of at least 1 − δPB − δΣ . The above bound is no longer an oracle bound but a data-driven one. A valid ϵΣ may be obtained from standard Gaussian sample covariance concentration bounds, and we provide an explicit choice in Appendix D.6. For a common choice of a diagonal Gaussian posterior, Σρ = diag(σ12 , . . . , σp2 ), Pp the mismatch term becomes ϵΣ k=1 σk2 ∥Mk ∥2F . Thus, the additional derandomization term is a weighted posterior-variance penalty, with weights determined entirely by the SLS basis directions. Moreover, even in the general case, MG can be computed offline.

6

Numerical experiments

All experiments are conducted on a laptop with a 12th Gen Intel i7-12700H processor, 32 GB of RAM, running Ubuntu 24.04.4 LTS. Code is implemented in Julia [Bezanson et al., 2017] (version 1.12.6) using automatic differentiation via ForwardDiff [Revels et al., 2016]. 6.1

Posterior learning versus point optimization in the low-data regime

We consider the double-integrator system as in (6) with     1 1 0.5 A= , B= , Q = diag(1, 0.1), 0 1 1

P = diag(2, 0.2),

R = 0.5,

with horizon T = 10 and Gaussian disturbance trajectory w = (x̄0 , w̄0 , . . . , w̄T −1 ) ∼ N (0, 0.1I). PT −1 PT −1 ⊤ ⊤ The finite-horizon quadratic cost is ℓ(x, u) = t=0 x̄⊤ t Qx̄t + t=0 ūt Rūt + x̄T P x̄T . The SLS parameterization has nullspace dimension p = 110. The loss in SLS coordinates is ℓ(α, w) = ∥M (α)w∥22 . Moreover, assume access to a dataset S = {wi }ni=1 . 10

Figure 1: Double-integrator low-data sweep with p = 110. Left: held-out test cost. Right: deterministic deployed sensitivity. Raw ERM is highly unstable when n ≪ p. Ridge regularization stabilizes point optimization, but the data-driven derandomized PAC-Bayes mean response achieves substantially lower test cost and lower sensitivity in the low-data regime. As n increases, all methods converge, indicating that the advantage is primarily finite-sample regularization. Low-data mean-response deployment. Next, we test whether the deterministic mean response obtained from a PAC-Bayes posterior acts as a useful finite-sample regularizer. We use the doubleintegrator system described above. The SLS nullspace has dimension p = 110. Thus, for small sample sizes n ≪ p, unregularized point optimization is expected to be statistically fragile. Hence, we treat the derandomized PAC-Bayes version as a synthesis method for control algorithms. We compare four deterministic deployments. The raw ERM baseline solves b S (α), αERM ∈ arg min L α

while ridge (Tikhonov) ERM solves b S (α) + γ∥α∥22 . αridge ∈ arg min L α

The proposed method optimizes a diagonal Gaussian posterior ρ = N (µρ , diag(σ12 , · · · , σp2 )) using the data-driven deterministic mean-response certificate b S (µρ ) + Compdata (ρ) + εΣ L

p X

σi2 ∥Mi ∥2F ,

i=1

and deploys the deterministic response µρ . As a prescient (clairvoyant) reference, we also report an oracle derandomized PAC-Bayes baseline that uses the true covariance Σw and population curvature G inside the certificate. This oracle is not available in a data-driven setting, but it indicates how much conservatism is introduced by covariance estimation. The prior is set to N (0, I). Moreover, for the epsilon term εΣ , we need a covariance scale, which is assumed to be σ̄ 2 = 0.1. All methods are evaluated on Ntest = 5000 newly drawn disturbance trajectories. We report the held-out test cost and the deterministic closed-loop sensitivity coefficient C(α) := ∥M (α)Σw M (α)⊤ ∥2F = ∥Σy (α)∥2F , the squared Frobenius norm of the induced closed-loop covariance, evaluated at the deployed response (the point solution for ERM/ridge, and the mean response µρ for the PAC-Bayes deployments). This is the quantity whose posterior average is the sensitivity coefficient Cρ = Eα∼ρ [C(α)] used in the Gaussian certificate. The advantage is most obvious when n < p. At n = 10 and n = 20, raw ERM often finds responses that incur very large held-out variance. Ridge regularization removes severe failures but remains substantially worse than the data-driven PAC-Bayes mean response. The data-driven derandomized method beats ridge ERM on every seed through n = 100, and the gap closes as n grows. This is the expected behavior. With enough trajectories, point optimization becomes well-conditioned, while in the low-data regime, the PAC-Bayes posterior acts as a sensitivity-aware regularizer whose mean response gives a stable deterministic controller. 11

Table 1: Double-integrator low-data sweep. Results are mean ± one standard deviation over 20 seeds. The PAC-Bayes columns report deterministic deployment of the posterior mean response. Held-out test cost n

raw ERM

ridge ERM

data PB derand oracle PB derand

10 20 40 100 400

164.01 ± 489.35 84.37 ± 216.30 4.08 ± 0.19 3.40 ± 0.09 3.19 ± 0.02

4.96 ± 0.33 4.47 ± 0.35 3.69 ± 0.11 3.32 ± 0.07 3.18 ± 0.02

3.22 ± 0.02 3.23 ± 0.02 3.24 ± 0.02 3.23 ± 0.03 3.18 ± 0.02

n

raw ERM

3.28 ± 0.03 3.27 ± 0.02 3.24 ± 0.03 3.20 ± 0.03 3.16 ± 0.02

Deployed sensitivity C(α)

10 2.42×105 ± 1.03×106 20 5.39×104 ± 2.07×105 40 2.11 ± 0.35 100 1.17 ± 0.11 400 0.97 ± 0.02

ridge ERM

data PB derand oracle PB derand

4.19 ± 1.01 3.20 ± 0.87 1.64 ± 0.18 1.14 ± 0.08 0.99 ± 0.02

0.85 ± 0.00 0.86 ± 0.01 0.87 ± 0.01 0.88 ± 0.01 0.88 ± 0.01

1.17 ± 0.13 1.11 ± 0.12 1.05 ± 0.10 0.97 ± 0.07 0.92 ± 0.01

Figure 2: Numerical certificate tightness for the learned data-driven posterior. Curves show means over 20 seeds. The prescient curve evaluates the same learned posterior with population quantities and is not used during training. The oracle derandomized baseline is useful as a reference but is not a deployable method since it uses population quantities that are unavailable from data. The fact that the data-driven method is close to this oracle in held-out cost, despite having a more conservative certificate, suggests that the learned response is very good and that the data-driven PAC-Bayes bound can be thought of as a good synthesis method for low-data learning control algorithms. The remaining gap is mainly a certification gap rather than a performance gap. Numerical certificates Beyond synthesis quality, we evaluate the numerical tightness of the datadriven deterministic mean-response certificate. For each sample size, we compare p X b Bdata (ρ) = LS (µρ ) + Compdata (ρ) + εΣ σk2 ∥Mk ∥2F k=1

with the held-out estimate of L(µρ ). We also report a prescient evaluation of the same learned posterior using the true disturbance covariance and population curvature. This diagnostic isolates the conservatism introduced by covariance estimation without changing the deployed mean response. The data-driven approach stays above the prescient one for all data points. Its average certified gap contracts from 94.29 at n = 10 to 2.49 at n = 400. At n = 10, the difference between the data-driven and prescient evaluations is 86.99, showing that covariance-estimation conservatism dominates in the lowest-data regime. At n = 400, this difference decreases to 1.32, while the prescient certificate remains 1.18 above the held-out cost. Thus, the certificate tightens steadily with data. 12

Table 2: Numerical tightness of the data-driven deterministic mean-response certificate. All quantities are evaluated at the posterior learned by the data-driven derandomized objective. Results are mean ± one standard deviation over 20 seeds. The prescient column evaluates the same learned posterior using the true covariance and population curvature; it is diagnostic and is not available to the learning algorithm. n held-out cost data certificate prescient certificate certificate − held-out 10 20 40 100 400

7

3.22 ± 0.02 3.23 ± 0.02 3.24 ± 0.02 3.23 ± 0.03 3.18 ± 0.02

97.51 ± 0.83 44.48 ± 0.83 22.57 ± 0.37 11.19 ± 0.19 5.68 ± 0.11

10.52 ± 0.20 8.23 ± 0.29 6.73 ± 0.16 5.52 ± 0.11 4.36 ± 0.07

94.29 ± 0.83 41.25 ± 0.82 19.34 ± 0.36 7.96 ± 0.19 2.49 ± 0.11

Conclusion

We developed PAC-Bayesian certificates for finite-horizon quadratic closed-loop control by placing posterior distributions directly over feasible SLS response coordinates. This parameterization preserves closed-loop achievability under unbounded posteriors and exposes the native loss as quadratic in both the disturbance trajectory and the free response coordinate. For Gaussian disturbances with arbitrary covariance this yields an exact one-sided Chernoff transform and a tractable sensitivity-based upper bound. The central message is that randomized certificates need not imply randomized deployment. The posterior serves as a finite-KL certificate-carrying object during learning, while the deployed controller is the deterministic mean response, which pays only the empirical-to-population curvature mismatch along the posterior covariance. Experiments on a double integrator confirm the resulting interpretation as a sensitivity-aware finite-sample regularizer, with lower held-out cost and sensitivity than ridge in the low-data regime, and the convergence of all methods as n grows. The sharpest certificates here rely on Gaussian or bounded disturbances and finite-horizon linear systems without hard state-input constraints. Extending the approach to constrained MPC and model uncertainty is left for future work.

References Pierre Alquier. User-friendly introduction to pac-bayes bounds. arXiv preprint arXiv:2110.11216, 2021. Anastasios N Angelopoulos, Rina Foygel Barber, and Stephen Bates. Theoretical foundations of conformal prediction. arXiv preprint arXiv:2411.11824, 2024. Michael Athans and Peter L Falb. Optimal control: an introduction to the theory and its applications. Courier Corporation, 2013. Alberto Bemporad, Manfred Morari, Vivek Dua, and Efstratios N Pistikopoulos. The explicit linear quadratic regulator for constrained systems. Automatica, 38(1):3–20, 2002. Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach to numerical computing. SIAM review, 59(1):65–98, 2017. URL https://doi.org/10.1137/ 141000671. Mahrokh Ghoddousi Boroujeni, Clara Lucia Galimberti, Andreas Krause, and Giancarlo FerrariTrecate. A pac-bayesian framework for optimal control with stability guarantees. In 2024 IEEE 63rd Conference on Decision and Control (CDC), pages 8237–8244, 2024. Mahrokh Ghoddousi Boroujeni, Clara Lucia Galimberti, Andreas Krause, and Giancarlo FerrariTrecate. Pac-bayesian optimal control with stability and generalization guarantees. arXiv preprint arXiv:2512.02858, 2025. Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press, Oxford, 2013. 13

Stephen P Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004. Marco C Campi and Simone Garatti. Introduction to the scenario approach. SIAM, 2018. Marco C Campi, Simone Garatti, and Maria Prandini. The scenario approach for systems and control design. Annual Reviews in Control, 33(2):149–157, 2009. Ioar Casado, Luis A Ortega, Aritz Pérez, and Andrés R Masegosa. Pac-bayes-chernoff bounds for unbounded losses. Advances in Neural Information Processing Systems, 37:24350–24374, 2024. Olivier Catoni. PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning. Institute of Mathematical Statistics, 2007. Kenneth R Davidson and Stanislaw J Szarek. Local operator theory, random matrices and banach spaces. In Handbook of the geometry of Banach spaces, volume 1, pages 317–366. Elsevier, 2001. Gintare Karolina Dziugaite and Daniel M. Roy. Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data. In Proceedings of the Thirty-Third Conference on Uncertainty in Artificial Intelligence, UAI 2017, Sydney, Australia, August 11-15, 2017. AUAI Press, 2017. Florian Dörfler, Jeremy Coulson, and Ivan Markovsky. Bridging direct and indirect data-driven control formulations via regularizations and relaxations. IEEE Transactions on Automatic Control, 68(2):883–897, 2023. doi: 10.1109/TAC.2022.3148374. Pascal Germain, Francis Bach, Alexandre Lacoste, and Simon Lacoste-Julien. Pac-bayesian theory meets bayesian inference. Advances in Neural Information Processing Systems, 29, 2016. Lars Grüne and Jürgen Pannek. Nonlinear model predictive control. In Nonlinear model predictive control: Theory and algorithms, pages 45–69. Springer, 2016. Benjamin Guedj. A primer on pac-bayesian learning. ArXiv, abs/1901.05353, 2019. David Lee Hanson and Farroll Tim Wright. A bound on tail probabilities for quadratic forms in independent random variables. The Annals of Mathematical Statistics, 42(3):1079–1083, 1971. Domagoj Herceg and Duarte Antunes. Distributionally robust pac-bayesian control. arXiv preprint arXiv:2604.10588, 2026. John Langford and Matthias Seeger. Bounds for averaging classifiers. 2001. Lars Lindemann, Yiqi Zhao, Xinyi Yu, George J Pappas, and Jyotirmoy V Deshmukh. Formal verification and control with conformal prediction. arXiv preprint arXiv:2409.00536, 2024. Anirudha Majumdar and Maxwell Goldstein. Pac-bayes control: Synthesizing controllers that provably generalize to novel environments. In Aude Billard, Anca Dragan, Jan Peters, and Jun Morimoto, editors, Proceedings of The 2nd Conference on Robot Learning, volume 87 of Proceedings of Machine Learning Research, pages 293–305. PMLR, 29–31 Oct 2018. Anirudha Majumdar, Alec Farid, and Anoopkumar Sonar. Pac-bayes control: Learning policies that provably generalize to novel environments, 2020. URL https://arxiv.org/abs/1806.04225. Ivan Markovsky and Florian Dörfler. Behavioral systems theory in data-driven analysis, signal processing, and control. Annual Reviews in Control, 52:42–64, 2021. doi: 10.1016/j.arcontrol. 2021.09.005. Arak M Mathai and Serge B Provost. Quadratic Forms in Random Variables: Theory and Applications. Marcel Dekker, Inc., New York, 1992. ISBN 0-8247-8691-2. Andreas Maurer. A note on the pac bayesian theorem. arXiv preprint cs/0411099, 2004. David A McAllester. Some pac-bayesian theorems. In Proceedings of the eleventh annual conference on Computational learning theory, pages 230–234, 1998. David A McAllester. Pac-bayesian stochastic model selection. Machine Learning, 51(1):5–21, 2003. 14

María Pérez-Ortiz, Omar Rivasplata, John Shawe-Taylor, and Csaba Szepesvári. Tighter risk certificates for neural networks. Journal of Machine Learning Research, 22(227):1–40, 2021. {Kaare Brandt} Petersen and {Michael Syskind} Pedersen. The Matrix Cookbook. Technical University of Denmark, 2006. Version 20051003. James B Rawlings, David Q Mayne, and Moritz M Diehl. Model predictive control: theory, computation, and design. (No Title), 2020. Benjamin Recht. A tour of reinforcement learning: The view from continuous control. Annual Review of Control, Robotics, and Autonomous Systems, 2:253–279, 2019. doi: 10.1146/ annurev-control-053018-023825. J. Revels, M. Lubin, and T. Papamarkou. Forward-mode automatic differentiation in Julia. arXiv:1607.07892 [cs.MS], 2016. URL https://arxiv.org/abs/1607.07892. Mark Rudelson and Roman Vershynin. Non-asymptotic theory of random matrices: extreme singular values. In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures, pages 1576–1602. World Scientific, 2010. Mark Rudelson and Roman Vershynin. Hanson-wright inequality and sub-gaussian concentration. 2013. René L Schilling. Measures, integrals and martingales. Cambridge University Press, 2017. Lukas Schüepp, Giulia De Pasquale, Florian Dörfler, and Carmen Amo Alonso. System level synthesis for affine control policies: Model-based and data-driven settings. In 2025 IEEE 64th Conference on Decision and Control (CDC), pages 1986–1992. IEEE, 2025. Jerome Sieber, Samir Bennani, and Melanie N Zeilinger. A system level approach to tube-based model predictive control. IEEE Control Systems Letters, 6:776–781, 2021. Leslie G Valiant. A theory of the learnable. Communications of the ACM, 27(11):1134–1142, 1984. P.C.N. Verheijen, V. Breschi, and M. Lazar. Handbook of linear data-driven predictive control: Theory, implementation and design. Annual Reviews in Control, 56:100914, 2023. ISSN 1367-5788. URL https://www.sciencedirect.com/science/article/pii/S1367578823000780. Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018. Eleftherios E Vlahakis, Lars Lindemann, Pantelis Sopasakis, and Dimos V Dimarogonas. Conformal prediction for distribution-free optimal control of linear stochastic systems. IEEE Control Systems Letters, 8:2835–2840, 2024. Martin J Wainwright. High-dimensional statistics: A non-asymptotic viewpoint, volume 48. Cambridge university press, 2019. Yuh-Shyang Wang, Nikolai Matni, and John C Doyle. A system-level approach to controller synthesis. IEEE Transactions on Automatic Control, 64(10):4079–4093, 2019. Ingvar Ziemann. An elementary proof of the hanson-wright inequality, 2025. URL https://arxiv. org/abs/2509.00881.

A

Extended Proofs

A.1

Proof of Proposition 5

Proof. Since yα = M (α)w + m(α) ∼ N (µy (α), Σy (α)), we can write

ℓ(α, w) = ∥yα ∥22 ,

L(α) = tr(Σy (α)) + ∥µy (α)∥22 . 15

Hence

i h i h 2 Ew eλ(L(α)−ℓ(α,w)) = eλL(α) Eyα e−λ∥yα ∥2 .

For a Gaussian vector y ∼ N (µ, Σ), with Σ ≻ 0 and any λ ≥ 0,   2 Ee−λ∥y∥2 = det(I + 2λΣ)−1/2 exp − λ µ⊤ (I + 2λΣ)−1 µ by [Mathai and Provost, 1992, Corollary 3.2a.2]. However, the same formula extends to case Σ ⪰ 0 almost trivially because of the negative sign in the exponent. The key observation is that I + 2λΣ is always well defined. Indeed, write the spectral decomposition of Σ ⪰ 0 as    Λ 0 Q⊤ R , Σ = [QR Q0 ] R 0 0 Q⊤ 0 where ΛR ≻ 0, QR spans Range(Σ), and Q0 spans Null(Σ). Equivalently, Σ = QR ΛR Q⊤ R, since the covariance has zero eigenvalues on Null(Σ). Writing ⊤ µ = QR Q⊤ R µ + Q0 Q0 µ = µR + µ0 ,

we have the representation   1/2 y = µ0 + QR Q⊤ µ + Λ z , R R

z ∼ N (0, Ir ).

Thus 1/2

∥y∥2 = ∥µ0 ∥2 + Q⊤ R µ + ΛR z

(30)

2

.

Applying the nonsingular formula on the r = rank(Σ) dimensional range space gives
2 −1 Ee−λ∥y∥ = det(Ir + 2λΛR )−1/2 exp −λ∥µ0 ∥2 − λµ̂⊤ µ̂R , R (Ir + 2λΛR ) where µ̂R = Q⊤ R µ is the vector in reduced basis. Since the eigenvalues of I + 2λΣ on the nullspace are equal to one, this is equivalent to the full-space expression
2 Ee−λ∥y∥ = det(I + 2λΣ)−1/2 exp −λµ⊤ (I + 2λΣ)−1 µ .

A.2

Proof of Corollary 6

Proof. We start from eq. (23) 1 Λα (λ) = λ tr(Σy (α))− log det(I +2λΣy (α))+λ∥µy (α)∥22 −λ µy (α)⊤ (I +2λΣy (α))−1 µy (α), 2 and group the terms with µy (α) and Σy (α) separately. For the mean term, we use [Petersen and Pedersen, 2006, Section 3.25, eq. (165)] which for our case reads as I − (I + 2λΣy (α))−1 = 2λΣy (α)(I + 2λΣy (α))−1 , to obtain λ∥µy (α)∥22 − λ µy (α)⊤ (I + 2λΣy (α))−1 µy (α) = 2λ2 µy (α)⊤ Σy (α)(I + 2λΣy (α))−1 µy (α) Note that

(I + 2λΣy (α))−1 ⪯ I because Σy (α) is PSD and λ ≥ 0. Hence 2λ2 µy (α)⊤ Σy (α)(I + 2λΣy (α))−1 µy (α) ≤ 2λ2 µy (α)⊤ Σy (α)µy (α). For the covariance term, let s1 , . . . , sr denote the eigenvalues of Σy (α). Then, due to basic properties of the trace and the determinant of a PSD matrix, we have  r  X 1 1 λ tr(Σy (α)) − log det(I + 2λΣy (α)) = λsi − log(1 + 2λsi ) . (31) 2 2 i=1 16

Next, we aim to find a simple bound for the above. It’s easy to recognize the first part of the scalar inequality, which we bound with a quadratic term, i.e. x−

1 log(1 + 2x) ≤ x2 , 2

∀x ≥ 0,

(32)

The validity of the bound can be checked by defining f (x) = x2 − x + 12 log(1 + 2x). We have 4x2 f (0) = 0, hence the inequality is tight at zero. Furthermore, f ′ (x) = 1+2x , which is clearly nonnegative for any x ≥ 0, proving our ad-hoc bound. It is also clear that the bound is tight around zero, but it gets progressively worse for higher values. However, optimization over λ keeps it in the "good" regime. Substituting x = λsi into (32) and comparing to (31), we obtain λ tr(Σy (α)) −

r X 1 log det(I + 2λΣy (α)) ≤ λ2 s2i = λ2 ∥Σy (α)∥2F , 2 i=1

(33)

by the definition of the Frobenius norm. Combining the two bounds proves (34). Define ψ(α, λ) := λ2 ∥Σy (α)∥2F + 2λ2 µy (α)⊤ Σy (α)µy (α)

(34)

Clearly ψ(α, 0) = ψ ′ (α, 0) = 0 as λ2 is the leading term in the expression, and moreover, it is convex. Finally, it is a valid object for Theorem 2 as stated in the paper which we inherited from [Casado et al., 2024].

B

Detailed affine SLS formulation

This appendix expands the finite-horizon affine System Level Synthesis Wang et al. [2019] (SLS) formulation used in the main text. The purpose is to clarify how the closed-loop response variables are defined and what achievability constraints they must satisfy. This is based on the known theory, and we mostly follow the exposition in [Schüepp et al., 2025] We restrict this appendix to the classical pointwise view of SLS, as the main purpose is to get readers unfamiliar with this framework up to speed. In addition, we provide vectorization of the constraints which are important for placing the posterior distribution over the response variables. In the subsequent Appendix we will show why the quadratic control loss becomes a quadratic function of the disturbance trajectory. B.1

Finite-horizon stacked dynamics

Consider the finite-horizon linear system x̄t+1 = Ax̄t + B ūt + w̄t ,

t = 0, . . . , T − 1,

where x̄t ∈ Rnx , ūt ∈ Rnu , and w̄t ∈ Rnx . We stack the state, input, and disturbance trajectories as     x̄0 ū0  x̄1   ū1  (T +1)nx T nu   x :=  , u :=  , (35)  ...  ∈ R  ...  ∈ R x̄T ūT −1 and  x̄0  w̄0  (T +1)nx  w :=  .  ...  ∈ R 

(36)

w̄T −1 Here w contains both the initial condition and the process disturbances. Thus, one sample wi in the PAC-Bayes formulation corresponds to one full finite-horizon disturbance trajectory as defined in (36). 17

Lifted system matrices.

For the stacked trajectories define   IT ⊗ B A := IT +1 ⊗ A, B := . 0nx ×T nu

Let Z ∈ R(T +1)nx ×(T +1)nx be the block down-shift operator. Then     0 0  Bu0   Ax0   Bu   Ax     1  1 ZAx =   , ZBu =  .  . .  ..   ..  BuT −1 AxT −1 It is useful to define ZA := ZA,

ZB := ZB,

Therefore, the stacked dynamics are x = ZA x + ZB u + w, or equivalently (I − ZA )x − ZB u = w. B.2

Affine closed-loop responses

Rather than parameterizing the controller directly, SLS parameterizes the closed-loop maps from disturbance trajectories to state and input trajectories. In the affine case, we write x = Φx w + ϕx ,

u = Φu w + ϕ u ,

or, equivalently,       x Φx ϕ = w+ x . u Φu ϕu Here and

Φx ∈ R(T +1)nx ×(T +1)nx , ϕx ∈ R(T +1)nx ,

Φu ∈ RT nu ×(T +1)nx , ϕu ∈ RT nu .

The matrices Φx and Φu are required to be block lower triangular, reflecting causality. The state and input at time t may depend only on the initial condition and disturbances up to the appropriate time, not on future disturbances. Substituting the affine response into the stacked dynamics gives (I − ZA )(Φx w + ϕx ) − ZB (Φu w + ϕu ) = w. Since this equality must hold for all disturbance trajectories w, the linear and affine parts must match separately. Therefore, the SLS achievability constraints are (I − ZA )Φx − ZB Φu = I, and (I − ZA )ϕx − ZB ϕu = 0. The first equation states that the linear response from disturbances to trajectories is compatible with the dynamics. The second equation states that the affine offset is itself dynamically consistent. More generally, if the stacked dynamics include a deterministic affine offset d, (I − ZA )x − ZB u = w + d, where d denotes the stacked deterministic offset, then the affine achievability constraint becomes (I − ZA )ϕx − ZB ϕu = d. The purely linear dynamics considered above correspond to d = 0. Combining the linear and affine components, the affine SLS constraint can be written as   Φx ϕ x [I − ZA −ZB ] = [I d] . Φu ϕ u For the nominal linear system without a deterministic offset, d = 0. 18

B.3

Controller recovery

The SLS response variables describe closed-loop behavior. Under the usual finite-horizon causality conditions, they also recover an implementable affine feedback controller. Since Φx is causal and has an invertible causal structure [Schüepp et al., 2025, Eq. 7a], we may solve w = Φ−1 x (x − ϕx ). Substituting this into u = Φ u w + ϕu gives u = Φu Φ−1 x (x − ϕx ) + ϕu . Therefore, the corresponding affine controller is u = Kx + k, where K := Φu Φ−1 x ,

k := ϕu − Φu Φ−1 x ϕx .

Thus, any achievable affine SLS response induces an implementable causal affine controller. Conversely, any causal affine controller applied to the finite-horizon linear system induces closed-loop responses (Φx , Φu , ϕx , ϕu ) satisfying the SLS achievability constraints. In this sense, the SLS equations parameterize exactly the achievable finite-horizon closed-loop responses. B.4

Weighted trajectory and quadratic loss

Let the finite-horizon quadratic control cost be defined by positive semidefinite state weights and positive definite input weights. Define the stacked weighting matrices Q̄ := blkdiag(Q0 , Q1 , . . . , QT ),

R̄ := blkdiag(R0 , R1 , . . . , RT −1 ).

For example, one may take Qt = Q for t < T , QT = P , and Rt = R. Define the weighted closed-loop trajectory  1/2  Q̄ x y(x, u) := . R̄1/2 u Substituting the affine SLS response gives  y(θ, w) =

 Q̄1/2 (Φx w + ϕx ) . R̄1/2 (Φu w + ϕu )

Therefore y(θ, w) = M (θ)w + m(θ), where M (θ) :=

 1/2  Q̄ Φx , R̄1/2 Φu

m(θ) :=

 1/2  Q̄ ϕx . R̄1/2 ϕu

The finite-horizon quadratic trajectory cost is then ℓ(θ, w) := ∥y(θ, w)∥22 = ∥M (θ)w + m(θ)∥22 . Expanding, ℓ(θ, w) = w⊤ M (θ)⊤ M (θ)w + 2m(θ)⊤ M (θ)w + m(θ)⊤ m(θ). Equivalently, defining Mc (θ) := M (θ)⊤ M (θ),

c(θ) := M (θ)⊤ m(θ),

we have ℓ(θ, w) = w⊤ Mc (θ)w + 2c(θ)⊤ w + m(θ)⊤ m(θ). This is the key structural property used in the main text: under affine SLS, the closed-loop quadratic control cost is a quadratic function of the disturbance trajectory w. 19

B.5

Vectorized SLS constraints

We vectorize the SLS response variables to expose the affine feasible set as a standard finitedimensional linear system Hϑ = h. This allows us to parameterize all feasible responses by a particular solution plus a nullspace basis, ϑ = ϑ0 + N α. Moreover, it allows us to place PACBayesian priors and posteriors directly on the free coordinates α. We now describe how the affine SLS achievability constraints can be written as a finite-dimensional linear system. Recall the affine closed-loop response x = Φx w + ϕx ,

u = Φu w + ϕ u ,

with achievability constraints (I − ZA )Φx − ZB Φu = I, and (I − ZA )ϕx − ZB ϕu = d, where d = 0 for the nominal linear system without an additional deterministic offset. Let F := [I − ZA

−ZB ] . Then the two achievability equations become     Φx ϕ F = I, F x = d. ϕu Φu

Define



 Φx Φ := , Φu



ϕ ϕ := x ϕu



to obtain FΦ = I, and F ϕ = d. Using the identity vec(ABC) = (C ⊤ ⊗ A) vec(B), the constraint FΦ = I can be vectorized as (see [Petersen and Pedersen, 2006]) (I ⊗ F) vec(Φ) = vec(I). The affine-offset constraint is already linear Fϕ = d. Thus, if we define the full vector of response variables   vec(Φ) ϑ := , ϕ the achievability constraints can be written compactly as Hϑ = h, where     I ⊗F 0 vec(I) H := , h := . 0 F d The response matrices Φx and Φu are block lower triangular due to causality. Rather than vectorizing all entries and imposing additional zero constraints, one may vectorize only the causal entries. Let S be a selection matrix that maps the reduced vector of causal entries ϑΦ to the full vectorization: vec(Φ) = SϑΦ . Then the decision variable becomes

  ϑ ϑ := Φ , ϕ

and the reduced SLS constraints are Hϑ = h, with H :=

 (I ⊗ F)S 0

 0 , F

 h :=

 vec(I) . d

This is the finite-dimensional linear system used for the nullspace parameterization. If ϑ0 is any solution of Hϑ = h and N is a basis for ker(H), then every feasible causal affine SLS response is represented as ϑ(α) = ϑ0 + N α. 20

(37)

Causality and selection matrices. The selection matrix S is only a helper matrix to remove entries of Φx and Φu that are structurally zero due to causality. Equivalently, one may vectorize all entries and add linear equality constraints enforcing the noncausal blocks to be zero. Both approaches lead to a linear system of the form Hϑ = h, and therefore both admit the same nullspace parameterization.

C

Mean-response deployment and quadratic loss structure

The goal of this section is to show that the affine nature of ϑ implies affine dependence in M (α), which induces quadratic L(α). Hence, the quadratic structure of L(α) is a consequence, not an assumption. Then we show how this ties into controller derandomization via mean deployment. We start in a somewhat backward order by first showing a more general statement that a quadratic L(α) admits a natural randomization. Proposition 9. Let L(α) = α⊤ Gα + 2g ⊤ α + α0 with G = G⊤ . Assume that ρ has a finite second moment and denote Eρ [α] = µρ , Cov(α) = Σρ . Then Eα∼ρ [L(α)] = L(µρ ) + tr(GΣρ ). Proof. Rewrite L(α) = tr(α⊤ Gα) + 2g ⊤ α + α0 . Due to the cyclicity of the trace, we have tr(α⊤ Gα) = tr(Gαα⊤ ). The linearity of both expectation and trace gives Eρ [tr(Gαα⊤ )] = tr(GEρ [αα⊤ ]). Next, we characterize Eρ [αα⊤ ]. By definition, Cov(α) := Eρ [(α − µρ )(α − µρ )⊤ ] = Eρ [αα⊤ ] − µρ µ⊤ ρ after simple algebraic manipulations. Moreover, it follows that Eρ [αα⊤ ] = Σρ + µρ µ⊤ . We have all the ingredients to compute Eρ [L(α)] = tr(GEρ [αα⊤ ]) + ρ ⊤ ⊤ ⊤ 2g ⊤ Eρ [α] + α0 = tr(G(Cov(α) + µρ µ⊤ ρ )) + 2g µρ + α0 = tr(GΣρ ) + tr(Gµρ µρ ) + 2g µρ + α0 . ⊤ ⊤ ⊤ Again, using the cyclicity of the trace, we have tr(Gµρ µρ ) = tr(µρ Gµρ ) = µρ Gµρ . The last equality is because µ⊤ ρ Gµρ is a scalar. Finally, recognizing the original definition of L(α) and plugging in µρ , we obtain the final form. Our next goal is to show that the SLS parametrization naturally results in a convex quadratic function in α, and the above proposition is readily applicable. C.1

From nullspace coordinates to affine cost maps

We make explicit how the nullspace parameterization of the SLS constraints induces affine dependence of the cost maps in the free coordinate α. Let dx := (T + 1)nx ,

du := T nu ,

Then

x ∈ Rdx , u ∈ Rdu , The SLS response variables have dimensions Φx ∈ Rdx ×dw ,

dw := (T + 1)nx . w ∈ Rdw .

Φu ∈ Rdu ×dw ,

and

ϕx ∈ Rdx , ϕu ∈ Rdu . We collect the response variables into the vector   vec(Φx ) vec(Φu ) ϑ :=  ∈ Rdϑ ϕx  ϕu where dϑ = dx dw + du dw + dx + du . The vectorized SLS constraints can be written as Hϑ = h. Let ϑ0 be one feasible solution, and let N = [N1

···

Np ] ∈ Rdϑ ×p

be a basis for ker(H). Then every feasible response can be written as ϑ(α) = ϑ0 + N α, 21

α ∈ Rp .

Equivalently, ϑ(α) = ϑ0 +

p X

αk Nk .

(38)

k=1

We partition the particular solution as   [0] vec(Φx )  [0]  vec(Φu ) ϑ0 =  . [0]  ϕx  [0] ϕu Similarly, each nullspace direction is partitioned as   [k] vec(Φx )  [k]  vec(Φu ) Nk =  ,  ϕ[k]  x [k] ϕu [k]

k = 1, . . . , p.

(39)

(40)

[k]

Here Φx and Φu are the full response matrices with their respective dimensions dx ×dw Φ[k] , x ∈R

du ×dw Φ[k] . u ∈R

Taking the corresponding blocks of ϑ(α) gives vec(Φx (α)) = vec(Φ[0] x )+ vec(Φu (α)) = vec(Φ[0] u )+

p X k=1 p X

αk vec(Φ[k] x ), αk vec(Φ[k] u ),

k=1

ϕx (α) = ϕ[0] x + ϕu (α) = ϕ[0] u +

p X k=1 p X

(41)

αk ϕ[k] x , αk ϕ[k] u .

k=1

Applying vec−1 to first two identities yields Φx (α) = Φ[0] x + Φu (α) = Φ[0] u +

p X k=1 p X

αk Φ[k] x ,

(42)

αk Φ[k] u .

(43)

k=1

Therefore, each SLS response component is affine in α. Now define the weighted cost maps  1/2   1/2  Q̄ Φx (α) Q̄ ϕx (α) M (α) := , m(α) := . R̄1/2 Φu (α) R̄1/2 ϕu (α) Let dy := dx + du . Then M (α) ∈ Rdy ×dw , m(α) ∈ Rdy . Substituting the affine response expansion gives " # " # p [0] [k] X Q̄1/2 Φx Q̄1/2 Φx M (α) = αk [0] + [k] , R̄1/2 Φu R̄1/2 Φu k=1 " # " # p [0] [k] X Q̄1/2 ϕx Q̄1/2 ϕx m(α) = αk [0] + [k] . R̄1/2 ϕu R̄1/2 ϕu k=1 22

(44)

(45)

Define # " [k] Q̄1/2 Φx Mk := [k] , R̄1/2 Φu

"

# [0] Q̄1/2 Φx M0 := [0] , R̄1/2 Φu

(46)

and "

# [0] Q̄1/2 ϕx m0 := [0] , R̄1/2 ϕu

"

# [k] Q̄1/2 ϕx mk := [k] , R̄1/2 ϕu

(47)

to obtain M (α) = M0 + m(α) = m0 +

p X k=1 p X

αk Mk ,

(48)

αk mk .

(49)

k=1

Thus M (α) and m(α) are affine functions of the SLS coordinate α. C.2

Quadratic loss in the SLS coordinate

We now show that the finite-horizon quadratic control loss is a convex quadratic function of α for each fixed disturbance trajectory w. For fixed w ∈ Rdw , define a0 (w) := M0 w + m0 ∈ Rdy , and

A(w) := [M1 w + m1 M2 w + m2 · · · Mp w + mp ] ∈ Rdy ×p . Using the affine expansions of M (α) and m(α), we obtain ! p p X X M (α)w + m(α) = M0 + αk Mk w + m0 + αk mk k=1

= M0 w + m0 + | {z } a0 (w)

(50)

k=1 p X

αk (Mk w + mk )

(51)

k=1

|

{z

A(w)α

}

= a0 (w) + A(w)α.

(52)

For a fixed w, we have ℓ(α, w) = ∥M (α)w + m(α)∥22 = ∥A(w)α + a0 (w)∥22 = α⊤ A(w)⊤ A(w)α + 2a0 (w)⊤ A(w)α + a0 (w)⊤ a0 (w). Hence, for every fixed w, the map α 7→ ℓ(α, w) is a quadratic function. Here we also note a biquadratic structure. The loss ℓ(w, α) is quadratic in α for fixed w and quadratic in w for a fixed α. Finally, assuming the relevant expectations are finite, the population risk is L(α) := Ew∼D [ℓ(α, w)] = α⊤ Gα α + 2gα⊤ α + cα , (53)     where Gα := Ew∼D  A(w)⊤ A(w) ∈ Rp×p , gα := Ew∼D A(w)⊤ a0 (w) ∈ Rp , and cα := Ew∼D a0 (w)⊤ a0 (w) ∈ R. Moreover, Gα ⪰ 0 ,thus L(α) is a convex quadratic function in the SLS coordinate α. C.3

Deterministic deployment (derandomization)

The PAC-Bayes certificate is naturally stated for the randomized posterior risk Eα∼ρ [L(α)]. In many control applications, it is highly preferable to deploy a single deterministic controller. This usually stems from certification and testing requirements. In the SLS quadratic setting, the natural choice is to deploy the posterior mean response. This subsection shows that the mean response is feasible, 23

inherits the PAC-Bayes certificate, and induces an implementable controller. This is a consequence of the convexity of the cost in response space, but we also show how it can be improved beyond the simple application of Jensen’s inequality due to quadratic structure. Jensen’s gap is exactly tr(Gα Σρ ). Let ρ be a posterior distribution over α. Assume that ρ has a finite second moment and define its mean and covariance by µρ := Eα∼ρ [α], Σρ := Covα∼ρ (α). Because ϑ(α) is affine in α, the posterior mean response is exactly the response obtained by evaluating the SLS parameterization at the posterior mean: ϑ̄ := Eα∼ρ [ϑ(α)] = ϑ(µρ ). Equivalently, denote Φ¯x = Φx (µρ ),

Φ¯u = Φu (µρ ),

ϕ¯x = ϕx (µρ ),

ϕ¯u = ϕu (µρ ).

Since the SLS achievability constraints are linear, feasibility is preserved under averaging. The key advantage of the quadratic SLS setting is that the population risk is convex quadratic in α as shown in the previous subsection. Therefore, Eα∼ρ [L(α)] = L(µρ ) + tr(Gα Σρ ). Since Gα ⪰ 0, the randomization gap (or tax) satisfies tr(Gα Σρ ) ≥ 0. Consequently, any PAC-Bayes upper bound on the randomized posterior risk also certifies the deterministic mean response. This can also be recovered from Jensen as a one-liner [Alquier, 2021] , but we want to explicitly characterize tr(GΣρ ) as a special important structure under quadratic costs. The term tr(GΣρ ) has a useful interpretation. It is the excess risk incurred by randomized deployment relative to mean-response deployment. It can also be seen as a curvature-weighted posterior variance term: posterior spread is expensive in directions where the closed-loop cost has high curvature and cheap in directions where the cost is flat. Why the Dirac-posterior issue does not arise? At first sight, certifying the deterministic mean response may appear to conflict with the usual PAC-Bayes requirement that the posterior have finite KL divergence with respect to the prior. Indeed, if the prior π is a continuous distribution, then the Dirac measure δµρ satisfies KL(δµρ ∥π) = ∞. Thus, one cannot usually apply the PAC-Bayes theorem directly to the deterministic posterior δµρ . The mean-response certificate does not do this. Instead, the PAC-Bayes theorem is applied to the non-degenerate posterior ρ, for which KL(ρ∥π) < ∞. Hence the deterministic mean response is certified indirectly through the randomized posterior certificate. The Dirac measure δµρ is never used as the PAC-Bayes posterior. In the quadratic SLS setting, this transfer is even more explicit. Since L(α) = α⊤ Gα + 2g ⊤ α + α0 ,

G ⪰ 0,

we have the exact identity Eα∼ρ [L(α)] = L(µρ ) + tr(GΣρ ). C.4

Recovering the controller for the mean response

The deterministic response certified by Proposition 7 is the mean SLS response (Φ¯x , Φ¯u , ϕ¯x , ϕ¯u ) = (Φx (µρ ), Φu (µρ ), ϕx (µρ ), ϕu (µρ )). To implement this response, we recover the corresponding affine controller using the standard SLS inversion. Since x = Φ¯x w + ϕ¯x , u = Φ¯u w + ϕ¯u , ¯ and Φx is causally invertible, we have −1 w = Φ¯x (x − ϕ¯x ).

Substituting this into the input response gives −1 u = Φ¯u Φ¯x (x − ϕ¯x ) + ϕ¯u .

24

Therefore, the deterministic affine controller that realizes the mean response is u = K̄x + k̄, where −1 −1 K̄ := Φ¯u Φ¯x , k̄ := ϕ¯u − Φ¯u Φ¯x ϕ¯x . Equivalently, K̄ = Φu (µρ )Φx (µρ )−1 , and k̄ = ϕu (µρ ) − Φu (µρ )Φx (µρ )−1 ϕx (µρ ). Thus, the deployment procedure is to compute µρ , then compute, (Φx (µρ ), Φu (µρ ), ϕx (µρ ), ϕu (µρ )) which finally give (K̄, k̄). This is not the same as averaging the feedback operators themselves. In general,   Eα∼ρ Φu (α)Φx (α)−1 ̸= Eα∼ρ [Φu (α)] Eα∼ρ [Φx (α)]−1 . The correct deterministic controller is obtained by first averaging in the response space and then recovering the controller.

D

Data-driven Gaussian certificates

This appendix shows how the Gaussian PAC-Bayes certificate can be made data-driven when the disturbance covariance Σw is unknown, hence transitioning from an oracle bound to an actually computable one. The key idea is to estimate Σw from the disturbance samples, inflate the estimate, and then use the inflated covariance in the Gaussian sensitivity coefficient. We also show that the same covariance event controls the curvature-mismatch term appearing in the deterministic mean response certificate. Finally, our approach maintains the curvature awareness interpretation. For simplicity, we focus on the zero offset setting (m(α) = 0) and zero mean disturbance, i.e. i.i.d.

wi ∼ N (0, Σw ), for every i ∈ {1, · · · , n}. The loss is ℓ(α, w) = ∥M (α)w∥2 , where M (α) = M0 +

p X

αk Mk .

k=1

Here α ∈ Rp , w ∈ Rdw , and Mk ∈ Rdy ×dw as explained in Appendix C.1 and is the consequence of the nullspace parametrization. D.1

Covariance event and inflated data-driven covariance

We can compute the empirical covariance as n X bw = 1 Σ wi wi⊤ n i=1

We assume that, for some ϵΣ > 0, the covariance event n o b w − Σw ∥op ≤ ϵΣ EΣ = ∥Σ

(54)

(55)

holds. We will discuss how to compute such a radius later on in Appendix D.6. Instead, we now focus on the implications of this event. Covariance event given in (55) implies b w + ϵΣ I. Σw ⪯ Σ The inflated empirical covariance is then b w + ϵΣ I. Σ Moreover, on the event EΣ , for every response α, it is not hard to see that b w + ϵΣ I)M (α)⊤ . M (α)Σw M (α)⊤ ⪯ M (α)(Σ

(56) (57)

Consequently b w + ϵΣ I)M (α)⊤ ∥2F . ∥M (α)Σw M (α)⊤ ∥2F ≤ ∥M (α)(Σ 25

(58)

D.2

Data-driven Gaussian complexity coefficient

The Gaussian PAC-Bayes certificate involves the oracle sensitivity coefficient, which we denote   Cρ = Eα∼ρ ∥M (α)Σw M (α)⊤ ∥2F . (59) Since Σw is unknown, Cρ is not directly computable (hence, oracle). On the covariance event EΣ , (58) gives the computable upper bound bρ , Cρ ≤ C where the data-driven sensitivity coefficient is h i bρ := Eα∼ρ ∥M (α)(Σ b w + ϵΣ I)M (α)⊤ ∥2 . C F

(60)

Thus, the unknown covariance is handled by replacing the oracle coefficient Cρ with the inflated bρ . For example, if the Gaussian PAC-Bayes theorem gives a bound of the form empirical coefficient C r Cρ (KL(ρ∥π) + log(n/δPB )) b Eα∼ρ [L(α)] ≤ Eα∼ρ [LS (α)] + 2 , n−1 then, on the covariance event, it yields the data-driven bound s b b S (α)] + 2 Cρ (KL(ρ∥π) + log(n/δPB )) . Eα∼ρ [L(α)] ≤ Eα∼ρ [L n−1 It is useful to define a data-driven complexity s Compdata (ρ) := 2

bρ (KL(ρ∥π) + log(n/δPB )) C n−1

and state the following. Proposition 10 (Data-driven randomized Gaussian certificate). Suppose the PAC-Bayes event holds bρ be as in with a probability of at least 1 − δPB and the covariance event EΣ holds. Let Cρ and C (59), (60). Then, on the intersection of these events, for all posteriors ρ ≪ π, b S (α)] + Compdata (ρ). Eα∼ρ [L(α)] ≤ Eα∼ρ [L

(61)

If P(EΣ ) ≥ 1 − δΣ , then the bound (61) holds with a probability of at least 1 − δPB − δΣ . Proof. The oracle Gaussian PAC-Bayes theorem gives, on the PAC-Bayes event, r b S (α)] + 2 Cρ (KL(ρ∥π) + log(n/δPB )) . Eα∼ρ [L(α)] ≤ Eα∼ρ [L n−1 bρ . Substitution gives the stated bound. The probability statement follows On EΣ , it holds that Cρ ≤ C by the union bound. D.3

Curvature matrices and mean-response decomposition

For fixed w, define the sample-dependent design matrix in the response coordinate α by A(w) := [M1 w

···

Mp w] ∈ Rdy ×p .

Then M (α)w = M0 w + A(w)α. Let a0 (w) := M0 w. The loss can be written as

ℓ(α, w) = ∥a0 (w) + A(w)α∥2 . 26

Expanding in α gives ℓ(α, w) = α⊤ A(w)⊤ A(w)α + 2a0 (w)⊤ A(w)α, +∥a0 (w)∥2 . From previous sections (see Appendix C.2), we know that     Gα := E A(w)⊤ A(w) , gα := E A(w)⊤ a0 (w) ,

  cα := E ∥a0 (w)∥2

and we can define the empirical counterparts as n X b S := 1 A(wi )⊤ A(wi ), G n i=1

Then

n

gbS :=

1X A(wi )⊤ a0 (wi ), n i=1

n

b cS :=

1X ∥a0 (wi )∥2 . n i=1

L(α) = α⊤ Gα α + 2gα⊤ α + cα ,

and

b S (α) = α⊤ G b S α + 2b L gS⊤ α + b cS .

Let

  Σρ := Eα∼ρ (α − µρ )(α − µρ )⊤ .

µρ := Eα∼ρ [α], In light of Proposition 9

Eα∼ρ [L(α)] = L(µρ ) + tr(GΣρ ), and b S (α)] = L b S (µρ ) + tr(G b S Σρ ). Eα∼ρ [L Combining these identities with Proposition 10 yields   b S (µρ ) + Compdata (ρ) + tr (G b S − G)Σρ . L(µρ ) ≤ L The remaining task is to somehow control the curvature mismatch   b S − G)Σρ . tr (G D.4

Curvature mismatch as a covariance-estimation error

We now show that the curvature mismatch can be controlled using the same covariance event EΣ . It will be useful to define
rρ (w) := tr A(w)Σρ A(w)⊤ . By the cyclicity of the trace,
rρ (w) = tr A(w)⊤ A(w)Σρ . Therefore,

n

b S Σρ ) = tr(G

1X rρ (wi ), n i=1

and tr(GΣρ ) = E[rρ (w)]. Hence

n   X b S − G)Σρ = 1 tr (G rρ (wi ) − E[rρ (w)]. n i=1

Next, we can expand rρ (w) since A(w) = [M1 w We have rρ (w) =

p X

···

Mp w] .

(Σρ )ij (Mi w)⊤ (Mj w).

i,j=1

27

Rewriting (Mi w)⊤ (Mj w) = w⊤ (Mi )⊤ Mj w, gives rρ (w) = w⊤ Aρ w, where Aρ :=

p X

(Σρ )ij (Mi )⊤ Mj ∈ Rdw ×dw .

(62)

i,j=1

Moreover Aρ ⪰ 0, because
2 w⊤ Aρ w = tr A(w)Σρ A(w)⊤ = ∥A(w)Σ1/2 ρ ∥F ≥ 0 for all w. Using w⊤ Aρ w = tr(Aρ ww⊤ ), we get n

1X b w ), rρ (wi ) = tr(Aρ Σ n i=1 and E[rρ (w)] = tr(Aρ Σw ). Therefore, we obtain the key identity     b S − G)Σρ = tr Aρ (Σ b w − Σw ) . tr (G

(63)

Lemma 11 (Curvature mismatch under covariance inflation). On the covariance event EΣ ,   b S − G)Σρ ≤ ϵΣ tr(Aρ ). tr (G b w − Σw ⪯ ϵΣ I. Since Aρ ⪰ 0, Proof. On EΣ , it holds Σ   b w − Σw ) ≤ tr (Aρ (ϵΣ I)) = ϵΣ tr(Aρ ). tr Aρ (Σ Using the identity (63) gives the result. For diagonal posterior covariance, Σρ = diag(σ12 , . . . , σp2 ), the matrix Aρ simplifies to Aρ =

p X

σk2 (Mk )⊤ Mk .

k=1

Consequently, tr(Aρ ) =

p X

σk2 ∥Mk ∥2F .

k=1

Thus, the mismatch penalty becomes the computable weighted posterior-variance penalty ϵΣ tr(Aρ ) = ϵΣ

p X

σk2 ∥Mk ∥2F .

k=1

D.5

Final data-driven deterministic certificate

Finally, we now combine the data-driven randomized PAC-Bayes certificate with the curvaturemismatch bound to round off the story via the following proposition. Proposition 12 (Data-driven deterministic certificate). Assume the zero mean, zero offset Gaussian i.i.d. setting wi ∼ N (0, Σw ) with m(α) = 0. Suppose the PAC-Bayes event holds with probability at least 1 − δPB and the covariance event EΣ holds with probability at least 1 − δΣ . Let ρ ≪ π be a posterior with Eρ [∥α∥2 ] < ∞, mean µρ and covariance Σρ . Then, with probability at least 1 − δPB − δΣ , b S (µρ ) + Compdata (ρ) + ϵΣ tr(Aρ ), L(µρ ) ≤ L 28

Pp where Aρ := i,j=1 (Σρ )ij Mi⊤ Mj as in (62). Equivalently, using tr(Aρ ) = tr(Σρ MG ) with the Gramian (MG )ij = tr(Mi⊤ Mj ), the mismatch term recovers the form of Proposition (8). For a diagonal posterior covariance Σρ = diag(σ12 , . . . , σp2 ), b S (µρ ) + Compdata (ρ) + ϵΣ L(µρ ) ≤ L

p X

σk2 ∥Mk ∥2F .

(64)

k=1

Proof. On the intersection of the PAC-Bayes event and EΣ , Proposition 10 gives b S (α)] + Compdata (ρ). Eα∼ρ [L(α)] ≤ Eα∼ρ [L

(65)

Since Eρ [∥α∥2 ] < ∞, the second-moment identities of Proposition 9 apply, and with Gα = bS , E[A(w)⊤ A(w)] and its empirical counterpart G Eα∼ρ [L(α)] = L(µρ ) + tr(Gα Σρ ),

b S (α)] = L b S (µρ ) + tr(G b S Σρ ). Eα∼ρ [L

(66)

Substituting both identities and rearranging yields
b S − Gα )Σρ . b S (µρ ) + Compdata (ρ) + tr (G L(µρ ) ≤ L

b S − Gα )Σρ = tr Aρ (Σ b w − Σw ) , Lemma 11 applies and gives By the key identity (63), tr (G
b S − Gα )Σρ ≤ ϵΣ tr(Aρ ) on EΣ . tr (G Pp CombiningP the two inequalities gives the first claim. For diagonal Σρ , Aρ = k=1 σk2 Mk⊤ Mk , hence p tr(Aρ ) = k=1 σk2 ∥Mk ∥2F , which yields the diagonal posterior formula. The probability statement follows by the union bound. D.6

Explicit Gaussian covariance radius

The preceding results are stated conditionally on a valid covariance event. We now give one explicit choice of ϵΣ for Gaussian disturbances. Assume i.i.d.

wi ∼ N (0, Σw ), Let

2 ∥Σw ∥op ≤ σ̄w .

n X bw = 1 wi wi⊤ . Σ n i=1

Fix δΣ ∈ (0, 1) and define r dw 2 log(2/δΣ ) an,δ := + . n n A standard Gaussian sample covariance concentration bound gives, with probability at least 1 − δΣ ,
b w − Σw ∥op ≤ ∥Σw ∥op 2an,δ + a2n,δ . ∥Σ r

2 Therefore, using the scale bound ∥Σw ∥op ≤ σ̄w , one may take
2 ϵΣ = σ̄w 2an,δ + a2n,δ .

This explicit choice is conservative and depends on the ambient disturbance dimension dw . Sharper covariance concentration inequalities, for example effective-rank bounds, can be substituted without changing the PSD-inflation and curvature-mismatch arguments above. The certificate only requires a valid event of the form b w − Σw ∥op ≤ ϵΣ . ∥Σ Next we prove that this radius is a sound choice. 29

Proof. Let r := rank(Σw ) ≤ dw . Write the reduced spectral decomposition Σw = QR ΛR Q⊤ R, where QR ∈ Rdw ×r has orthonormal columns and ΛR ∈ Rr×r is positive definite. Define 1/2

S := QR ΛR ∈ Rdw ×r , so that Σw = SS ⊤ . Since wi ∼ N (0, Σw ), we may write zi ∼ N (0, Ir ).

wi = Szi , Let Z ∈ R

n×r

be the matrix whose i-th row is zi⊤ . Then   n X 1 ⊤ bw = 1 Σ wi wi⊤ = S Z Z S⊤, n i=1 n

and hence

 b w − Σw = S Σ

 1 ⊤ Z Z − Ir S ⊤ . n

Therefore, b w − Σw ∥op ≤ ∥S∥2 ∥Σ op

1 ⊤ Z Z − Ir n op

= ∥Σw ∥op

1 ⊤ Z Z − Ir . n op

Let sj (Z) be the j-th singular value of the matrix Z, and let the minimum and maximum values be denoted by smin and smax . By the Gaussian extreme singular-value inequality [Davidson and Szarek, 2 2001][Rudelson and Vershynin, 2010, Theorem 2.6,eq (2.3)], with probability at least 1 − 2 exp (− t2 ) for t ≥ 0 √ √ smax (Z) ≤ n + r + t (67) and smin (Z) ≥ Setting t =

p

√

n−

√

r − t.

2 log(2/δΣ ) gives the probability of 1 − δΣ and p √ √ smax (Z) ≤ n + r + 2 log(2/δΣ )

(68)

(69)

and

p √ √ smin (Z) ≥ n − r − 2 log(2/δΣ ). √ Divide both inequalities by n and define r r r 2 log(2/δΣ ) a(n, δΣ , r) = + , n n we have  1 ⊤ Z Z − Ir ≤ max (1 + a(n, δΣ , r))2 − 1, 1 − [1 − a(n, δΣ , r)]2+ . n op

(70)

The clamping to one is necessary because if a > 1 the Z is rank deficient and has zero singular values, but due to the identity matrix, the overall matrix must have a singular value of one. Since (1 + a)2 − 1 = 2a + a2 and

1 − [1 − a]2+ ≤ 2a + a2 ,

we combine both into 1 ⊤ Z Z − Ir ≤ 2a(n, δΣ , r) + (a(n, δΣ , r))2 . n op Finally, a safe (and conservative) choice for r is r = dw , so a(n, δΣ , r) ≤ a(n, δΣ , dw ), and 2 ∥Σw ∥op ≤ σ̄w . This proves the claim. 30

D.7

Remark on affine offsets

The clean curvature-mismatch formula above assumes m(α) = 0. If affine offsets are present, then   w M (α)w + m(α) = M̄ (α)w̄, w̄ := , 1 where M̄ (α) = [M (α) m(α)] . The same algebra can be repeated with the augmented second-moment matrix   Σw 0 ⊤ E[w̄w̄ ] = 0 1 in the zero-mean case. However, if the mean is unknown or the offset terms are estimated from data, then an additional concentration event for the mean or the augmented second moment is needed. For this reason, the main data-driven curvature-mismatch refinement is stated in the zero-mean, zero-offset setting.

31