Average Gramian Matrix in Control & Randomness

• Average Gramian Matrix is a structured approach that aggregates trajectory-based inner products to capture controllability, observability, and energy functionals in dynamic systems.
• It leverages time, perturbation, ensemble, and parameter averaging to extend classical LTI Gramians to nonlinear and time-varying systems.
• In randomized algorithms and spherical designs, it provides precise moment-matching and approximation guarantees useful for numerical linear algebra and system identification.

Searching arXiv for the cited papers and related terminology. “Average Gramian Matrix” denotes, most explicitly in the empirical Gramian framework, a Gramian assembled by averaging trajectory-based inner products over multiple axes of variability—time, perturbations, parameters, operating points, and ensembles. Closely related literatures study averaged resolvents of random Gram matrices, expectations of random Gram products, Monte Carlo averages of weighted outer products, and Gramian constraints that reproduce spherical moments. Across these settings, the common structure is a positive-semidefinite matrix built from pairwise products and then averaged in a manner dictated by dynamics, randomness, sampling, or symmetry (Himpe, 2016, Alt et al., 2016, Auguin et al., 2017, Holodnak et al., 2013, Martin et al., 2013, Waldron, 5 Nov 2025).

1. Terminological scope and underlying matrix objects

In linear system theory, a Gramian is an operator encoding controllability, observability, or input-output coupling. For an asymptotically stable LTI system

$\dot x(t)=Ax(t)+Bu(t),\qquad y(t)=Cx(t),$

the classic Gramians are the controllability Gramian $W_c$, observability Gramian $W_o$, and, for square systems, the cross Gramian $W_x$. They admit both operator equations and integral representations: $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$

$$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$

$$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$

For state-space symmetric systems with $A=A^T$ and $B=C^T$, $W_c$, $W_c$0, and $W_c$1 coincide, which is the structural basis for cross-Gramian-based balancing (Himpe, 2016).

In matrix analysis and probability, the relevant object is often a Gram matrix or Gram product. For a matrix $W_c$2, the Gram matrix is $W_c$3 or $W_c$4; for sampled vectors $W_c$5, the Gramian is the pairwise inner-product matrix $W_c$6; for random matrices $W_c$7, the Gram matrix $W_c$8 is positive semidefinite and its spectrum is analyzed through the resolvent $W_c$9. These are not interchangeable definitions, but they are connected by the same bilinear-product structure. A plausible implication is that the phrase “average Gramian matrix” should be read contextually: in control it refers to averaged trajectory covariances, in random-matrix theory to averaged spectral observables or expectations, and in geometric design theory to moment-matching Gramian statistics.

2. Classical system Gramians and their empirical extensions

The empirical Gramian framework extends system Gramians from asymptotically stable LTI systems to nonlinear, parametric, and time-varying systems. For a nonlinear parametric system

$W_o$0

empirical Gramians are computed from simulated trajectories by averaging inner products of centered state or output trajectories across perturbations. The perturbation sets are given by canonical input directions $W_o$1, input scales $W_o$2, initial-state directions $W_o$3, and initial-state scales $W_o$4; signs $W_o$5 and scale subdivisions are controlled by flags, while rotations are limited to $W_o$6. Centered trajectories are

$W_o$7

where the centering choice can be steady-state, final state, arithmetic time average, RMS, or mid-range (Himpe, 2016).

The empirical controllability Gramian is assembled from state responses to input perturbations,

$W_o$8

The empirical observability Gramian is assembled from output trajectories caused by initial-state perturbations,

$W_o$9

and the empirical cross Gramian for square systems averages both input and initial-state perturbations,

$W_x$0

$W_x$1

Parameter-space analogs include the sensitivity Gramian $W_x$2, the identifiability Gramian

$W_x$3

and, for square systems, the joint Gramian with cross-identifiability

$W_x$4

These constructions formalize how local perturbation responses are converted into a single matrix summarizing input-output coherence and identifiability over an operating region (Himpe, 2016).

3. “Average” in the empirical Gramian framework

In emgr, an “Average Gramian Matrix” is any Gramian assembled by averaging trajectory-based inner products over multiple axes of variability: time, perturbations, parameters, operating points, and ensembles. Time averaging is implemented through integrals $W_x$5, with discrete approximation by sums $W_x$6. Ensemble averaging appears explicitly in the prefactors $W_x$7, $W_x$8, and in the scale normalizations $W_x$9, $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$0, or $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$1. Parametric averaging is performed by supplying multiple parameter columns $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$2; emgr computes a Gramian for each $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$3 and returns

$$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$4

Operating points can be varied through $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$5 and $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$6, and weighted averaging can be applied in post-processing via

$$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$7

(Himpe, 2016).

This averaging is not merely notational. Centering reduces bias from offsets and makes empirical Gramians robust to nonzero steady-states and non-impulse excitation. For asymptotically stable systems and impulse input, empirical Gramians converge to LTI Gramians as $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$8; for unstable or time-varying systems, one instead uses finite $$A W_c + W_c A^T + B B^T = 0,\qquad W_c = \int_0^\infty e^{At}BB^Te^{A^Tt}\,dt,$$9 and interprets the result as a time-limited Gramian. The framework therefore treats “average” as a controlled aggregation protocol rather than as a single analytic formula.

The same construction also induces scalar summaries. Average controllability for LTI networks is commonly quantified by $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$0; with empirical $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$1, the same scalar metric extends to nonlinear and parametric systems over the chosen operating region, for example $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$2 or $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$3. Other summaries listed in the framework are $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$4, smallest eigenvalues $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$5, condition numbers $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$6, sums of neglected singular values, and diagonal entries measuring state input-output importance or coherence. This suggests that an average Gramian matrix often functions as a finite-dimensional surrogate for a family of operating-condition-dependent energy functionals.

4. Averaged resolvents and deterministic equivalents for random Gram matrices

For random-matrix theory, the central object is a random $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$7 matrix $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$8 with independent, centered entries and arbitrary variance profile $$A^T W_o + W_o A + C^T C = 0,\qquad W_o = \int_0^\infty e^{A^Tt}C^TCe^{At}\,dt,$$9, where $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$0. The Gram matrix is $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$1, and the averaged spectral observable is the empirical Stieltjes transform

$$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$2

The limiting deterministic quantity is the average

$$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$3

where $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$4 is the unique holomorphic solution with $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$5 of the Gram Dyson equation

$$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$6

Its average is the Stieltjes transform of a probability measure $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$7,

$$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$8

which generalizes the Marchenko–Pastur law for arbitrary variance profiles (Alt et al., 2016).

Under flatness of $$A W_x + W_x A = -BC,\qquad W_x = \int_0^\infty e^{At}BC\,e^{At}\,dt.$$9, a primitivity-type lower regularity condition, bounded moments, and an aspect ratio bounded away from $A=A^T$0 and $A=A^T$1, optimal entrywise and averaged local laws hold down to spectral scales $A=A^T$2. In the bulk, the averaged local law yields

$A=A^T$3

with overwhelming probability uniformly in the spectral domain $A=A^T$4; away from the support of $A=A^T$5,

$A=A^T$6

These are sharper than the entrywise bounds, which are of order $A=A^T$7 in the bulk and $A=A^T$8 away from support. The paper also proves bulk rigidity,

$A=A^T$9

and absence of eigenvalues away from $B=C^T$0 with optimal accuracy.

The edge behavior splits into a hard-edge square regime and a soft-edge properly rectangular regime. For $B=C^T$1 under block fully indecomposable variance assumptions, $B=C^T$2 has no atom at $B=C^T$3 and its density satisfies $B=C^T$4, so the density behaves like $B=C^T$5 at the hard edge. For $B=C^T$6 with lower bounded variances, there is a macroscopic gap between $B=C^T$7 and the lower edge of $B=C^T$8, and if $B=C^T$9 then $W_c$0 and $W_c$1 a.w.o.p. A major consequence is the deterministic-equivalent principle: $W_c$2 with the example

$W_c$3

as in the MIMO channel-capacity motivation. In this literature, the “average” attached to a Gram matrix is the trace-normalized resolvent and the spectral measure it converges to (Alt et al., 2016).

5. Expected Gram matrices and average characteristic data

A second meaning of average Gramian is literal expectation. For a $W_c$4 complex central Gaussian matrix $W_c$5 with independent entries $W_c$6, arbitrary positive variance profile, and Gram matrix $W_c$7, one has

$W_c$8

Hence

$W_c$9

and $W_c$00 because the entries are independent, zero-mean, and circularly symmetric. The trace and determinant moments are also explicit: $W_c$01

$W_c$02

Here the average Gramian matrix is exactly the expectation of the random Gram product, and it records expected per-receive-antenna signal power under an arbitrary variance profile (Auguin et al., 2017).

A more general expectation theory is developed for random Gram matrices generated from sampled columns. Let $W_c$03 or $W_c$04 be formed from $W_c$05 independent columns $W_c$06, each distributed as a random vector $W_c$07, and let

$W_c$08

If

$W_c$09

is the second-moment matrix and $W_c$10, then the expected determinant and expected permanent of $W_c$11,

$W_c$12

have exponential generating functions

$W_c$13

$W_c$14

The corresponding recursions are

$W_c$15

with $W_c$16. Since $W_c$17, $W_c$18 for $W_c$19, hence $W_c$20 in that regime. The expected coefficients of the characteristic and permanental polynomials factor through these sequences via

$W_c$21

This framework shows that, for independent sampled columns, average Gramian information can be compressed to the second-moment matrix $W_c$22 and its power traces, without assuming Gaussianity or independence of coordinates within a column (Martin et al., 2013).

6. Monte Carlo averages of outer products and approximation guarantees

In randomized numerical linear algebra, the average Gramian matrix is represented by a finite Monte Carlo average of weighted rank-one outer products. For $W_c$23 with columns $W_c$24, one samples indices $W_c$25 independently with replacement according to probabilities $W_c$26 and forms

$W_c$27

The estimator is unbiased: $W_c$28 This gives a literal average-over-samples approximation to the Gram matrix, with weights chosen to compensate for nonuniform sampling (Holodnak et al., 2013).

The same paper gives exact and probabilistic characterizations. If $W_c$29 has rank $W_c$30, selected columns $W_c$31, and diagonal weights $W_c$32, then

$W_c$33

equivalently $W_c$34 has orthonormal rows. In the special case $W_c$35, the unique diagonal weights are

$W_c$36

Thus exact average representation depends on the right singular vector geometry rather than solely on column norms.

For approximation, the error bounds are dimension-free in the sense that they depend on the stable rank

$W_c$37

or on $W_c$38, not on ambient matrix dimensions. With nearly-optimal probabilities $W_c$39, $W_c$40, one obtains

$W_c$41

with probability at least $W_c$42, provided

$W_c$43

An intrinsic-dimension variant replaces $W_c$44 by $W_c$45, and leverage-score sampling yields a rank-based bound

$W_c$46

The practical guidance given is to prefer norm-squared probabilities

$W_c$47

which minimize $W_c$48 and empirically yield smaller spectral errors than leverage scores when $W_c$49 is small (Holodnak et al., 2013).

7. Gramian averages in spherical designs and moment matching

For real and complex spherical designs, the Gramian is the matrix of pairwise inner products,

$W_c$50

and with weights $W_c$51 one may pass to the weighted Gram

$W_c$52

If a weighted spherical design is defined as a cubature rule for a unitarily invariant polynomial space $W_c$53, then any unitary image is also such a design, and spherical designs for $W_c$54 are determined up to unitary equivalence by their Gramian. The reproducing kernel has the form

$W_c$55

and the associated potential

$W_c$56

vanishes exactly for weighted spherical designs. In this setting, Gramian-based energies encode averaged polynomial moments of pairwise inner products (Waldron, 5 Nov 2025).

The paper makes the “average Gramian” viewpoint explicit through moment identities. For balanced sets,

$W_c$57

and for the second moment,

$W_c$58

Equivalently, the average outer product is isotropic. For unit weights $W_c$59 and $W_c$60, this is the tight-frame identity

$W_c$61

In Gramian terms, $W_c$62 has rank $W_c$63 with $W_c$64 nonzero eigenvalues all equal to $W_c$65, and the frame potential

$W_c$66

satisfies

$W_c$67

with equality for tight frames or designs.

The potentials used to characterize designs are themselves averages of Gramian entries under suitable nonlinear functions. Examples include

$W_c$68

for real half-designs,

$W_c$69

for complex $W_c$70-designs, and

$W_c$71

for complex projective $W_c$72-designs. The design condition is that the sample Gramian reproduces spherical averages such as

$W_c$73

This suggests a precise geometric meaning of average Gramian matrix: a Gramian whose entrywise statistics and spectral profile match the spherical averages dictated by symmetry, with isotropy at second order and kernel-based moment matching at higher orders (Waldron, 5 Nov 2025).