---
title: Average Gramian Matrix in Control & Randomness
url: https://www.emergentmind.com/topics/average-gramian-matrix
type: topic
---

# Average Gramian Matrix in Control & Randomness

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 [1611.00675, 1606.07353, 1705.05214, 1310.1502, 1309.2599, 2511.07452].

## 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_o\), and \(W_x\) coincide, which is the structural basis for cross-Gramian-based balancing [1611.00675].

In matrix analysis and probability, the relevant object is often a Gram matrix or Gram product. For a matrix \(A\), the Gram matrix is \(AA^T\) or \(A^*A\); for sampled vectors \(x_i\), the Gramian is the pairwise inner-product matrix \(\Gamma=[\langle x_i,x_j\rangle]\); for random matrices \(X\), the Gram matrix \(XX^*\) is positive semidefinite and its spectrum is analyzed through the resolvent \((XX^*-\zeta I)^{-1}\). 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
\[
\dot x = f(t,x,u,p),\qquad y=g(t,x,u,p),
\]
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 \(E_u=\{e^m\}\), input scales \(S_u=\{c_k\}\), initial-state directions \(E_x=\{\varepsilon^j\}\), and initial-state scales \(S_x=\{d_l\}\); signs \(\pm\) and scale subdivisions are controlled by flags, while rotations are limited to \(\pm I\). Centered trajectories are
\[
\Delta x(t)=x(t)-\bar x,\qquad \Delta y(t)=y(t)-\bar y,
\]
where the centering choice can be steady-state, final state, arithmetic time average, RMS, or mid-range [1611.00675].

The empirical controllability Gramian is assembled from state responses to input perturbations,
\[
\widehat W_C
=
\frac{1}{|S_u|}
\sum_{k=1}^{|S_u|}
\sum_{m=1}^{M}
\frac{1}{c_k^2}
\int_0^T \Psi^{km}(t)\,dt,
\qquad
\Psi^{km}(t)
=
(x^{km}(t)-\bar x^{km})(x^{km}(t)-\bar x^{km})^T.
\]
The empirical observability Gramian is assembled from output trajectories caused by initial-state perturbations,
\[
\widehat W_O
=
\frac{1}{|S_x|}
\sum_{l=1}^{|S_x|}
\frac{1}{d_l^2}
\int_0^T \Psi^l(t)\,dt,
\qquad
(\Psi^l)_{ij}(t)
=
(y^{li}(t)-\bar y^{li})^T(y^{lj}(t)-\bar y^{lj}),
\]
and the empirical cross Gramian for square systems averages both input and initial-state perturbations,
\[
\widehat W_X
=
\frac{1}{|S_u||S_x|M}
\sum_{k=1}^{|S_u|}
\sum_{l=1}^{|S_x|}
\sum_{m=1}^{M}
\frac{1}{c_k d_l}
\int_0^T \Psi^{klm}(t)\,dt,
\]
\[
(\Psi^{klm})_{ij}(t)
=
(x_i^{km}(t)-\bar x_i^{km})(y_m^{lj}(t)-\bar y_m^{lj}).
\]
Parameter-space analogs include the sensitivity Gramian \(W_{S,ii}:=\operatorname{tr}(W_{C,i})\), the identifiability Gramian
\[
W_I = W_P - W_M^T W_O^{-1} W_M,
\]
and, for square systems, the joint Gramian with cross-identifiability
\[
W_{\ddot I}= -\frac{1}{2}W_m^T(W_X+W_X^T)^{-1}W_m.
\]
These constructions formalize how local perturbation responses are converted into a single matrix summarizing input-output coherence and identifiability over an operating region [1611.00675].

## 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 \(\int_0^T(\cdot)\,dt\), with discrete approximation by sums \(h\sum \phi(t_n)\). Ensemble averaging appears explicitly in the prefactors \(1/|S_u|\), \(1/|S_x|\), and in the scale normalizations \(1/c_k^2\), \(1/d_l^2\), or \(1/(c_k d_l)\). Parametric averaging is performed by supplying multiple parameter columns \(pr\in\mathbb R^{P\times K}\); emgr computes a Gramian for each \(\theta^{(k)}\) and returns
\[
\bar W = \frac{1}{K}\sum_{k=1}^K W^{(k)}.
\]
Operating points can be varied through \(x_s\) and \(u_s\), and weighted averaging can be applied in post-processing via
\[
\bar W=\sum_{k=1}^K w_k W^{(k)},\qquad \sum w_k=1
\]
[1611.00675].

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 \(T\to\infty\); for unstable or time-varying systems, one instead uses finite \(T\) 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 \(\operatorname{tr}(W_c)\); with empirical \(W_c\), the same scalar metric extends to nonlinear and parametric systems over the chosen operating region, for example \(\operatorname{AC}\approx \operatorname{tr}(\widehat W_C)\) or \(\operatorname{tr}(\widehat W_C)/n\). Other summaries listed in the framework are \(\operatorname{tr}(\widehat W_X)\), smallest eigenvalues \(\lambda_{\min}\), condition numbers \(\kappa(W)\), 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 \(p\times n\) matrix \(X=(x_{ik})\) with independent, centered entries and arbitrary variance profile \(S=(s_{ik})\), where \(s_{ik}=\mathbb E x_{ik}^2\). The Gram matrix is \(XX^*\), and the averaged spectral observable is the empirical Stieltjes transform
\[
m_N(\zeta)=\frac{1}{p}\operatorname{Tr}(XX^*-\zeta I_p)^{-1}.
\]
The limiting deterministic quantity is the average
\[
\{m(\zeta)\}=\frac{1}{p}\sum_i m_i(\zeta),
\]
where \(m(\zeta)\in\mathbb C^p\) is the unique holomorphic solution with \(\Im m(\zeta)>0\) of the Gram Dyson equation
\[
-\frac{1}{m(\zeta)}=\zeta - S\frac{1}{1+S^Tm(\zeta)}.
\]
Its average is the Stieltjes transform of a probability measure \(\nu\),
\[
\{m(\zeta)\}=\int_{\mathbb R}\frac{1}{\omega-\zeta}\,\nu(d\omega),
\]
which generalizes the Marchenko–Pastur law for arbitrary variance profiles [1606.07353].

Under flatness of \(S\), a primitivity-type lower regularity condition, bounded moments, and an aspect ratio bounded away from \(0\) and \(\infty\), optimal entrywise and averaged local laws hold down to spectral scales \(\eta=\Im\zeta\gtrsim p^{-1+\varepsilon}\). In the bulk, the averaged local law yields
\[
|m_N(\zeta)-\{m(\zeta)\}|
\le
\frac{p^\varepsilon}{p\,\Im\zeta}
\]
with overwhelming probability uniformly in the spectral domain \(\mathbb D_\delta\); away from the support of \(\nu\),
\[
|m_N(\zeta)-\{m(\zeta)\}|
\le
\frac{p^\varepsilon}{p}.
\]
These are sharper than the entrywise bounds, which are of order \((p\,\Im\zeta)^{-1/2}\) in the bulk and \(p^{-1/2}\) away from support. The paper also proves bulk rigidity,
\[
\mathbb P\!\left(|\lambda_{i(\tau)}-\tau|\ge p^\varepsilon/p\right)\le C_{\varepsilon,D}p^{-D},
\]
and absence of eigenvalues away from \(\operatorname{supp}\nu\) with optimal accuracy.

The edge behavior splits into a hard-edge square regime and a soft-edge properly rectangular regime. For \(p=n\) under block fully indecomposable variance assumptions, \(\nu\) has no atom at \(0\) and its density satisfies \(\lim_{\omega\downarrow 0}\rho(\omega)\sqrt{\omega}\in(0,\infty)\), so the density behaves like \(\omega^{-1/2}\) at the hard edge. For \(|p/n-1|\ge d_*>0\) with lower bounded variances, there is a macroscopic gap between \(0\) and the lower edge of \(\operatorname{supp}\rho\), and if \(p>n\) then \(\vartheta=1-n/p\) and \(\dim\ker(XX^*)=p-n\) a.w.o.p. A major consequence is the deterministic-equivalent principle:
\[
\frac{1}{p}\operatorname{Tr}f(XX^*) \approx \int f(\omega)\,\nu(d\omega),
\]
with the example
\[
\frac{1}{p}\log\det(I+\sigma^{-2}XX^*)
\approx
\int \log(1+\sigma^{-2}\omega)\,\nu(d\omega),
\]
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 [1606.07353].

## 5. Expected Gram matrices and average characteristic data

A second meaning of average Gramian is literal expectation. For a \(2\times2\) complex central Gaussian matrix \(H\) with independent entries \(H_{ij}\sim \mathcal{CN}(0,\phi_{ij})\), arbitrary positive variance profile, and Gram matrix \(W=HH^\dagger\), one has
\[
W_{11}=|H_{11}|^2+|H_{12}|^2,\qquad
W_{22}=|H_{21}|^2+|H_{22}|^2,\qquad
W_{12}=H_{11}H_{21}^*+H_{12}H_{22}^*.
\]
Hence
\[
\mathbb E[W]=\operatorname{diag}(\phi_{11}+\phi_{12},\,\phi_{21}+\phi_{22}),
\]
and \(\mathbb E[W_{12}]=0\) because the entries are independent, zero-mean, and circularly symmetric. The trace and determinant moments are also explicit:
\[
\mathbb E[\operatorname{tr}(W)] = \sum_{i,j}\phi_{ij},\qquad
\operatorname{Var}(\operatorname{tr}(W)) = \sum_{i,j}\phi_{ij}^2,
\]
\[
M_{\operatorname{tr}}(t)=\prod_{i,j}(1-\phi_{ij}t)^{-1},
\qquad
\mathbb E[\det(W)]=\phi_{11}\phi_{22}+\phi_{12}\phi_{21}.
\]
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 [1705.05214].

A more general expectation theory is developed for random Gram matrices generated from sampled columns. Let \(A\in\mathbb R^{t\times n}\) or \(\mathbb C^{t\times n}\) be formed from \(n\) independent columns \(a_1,\dots,a_n\), each distributed as a random vector \(w\), and let
\[
G_n=A^TA \quad \text{or}\quad G_n=A^*A.
\]
If
\[
M_{ij}=\mathbb E(w_i w_j)
\]
is the second-moment matrix and \(t_k=\operatorname{tr}(M^k)\), then the expected determinant and expected permanent of \(G_n\),
\[
a_n=\mathbb E[\det(G_n)],\qquad p_n=\mathbb E[\operatorname{per}(G_n)],
\]
have exponential generating functions
\[
\sum_{n=0}^\infty a_n\frac{x^n}{n!}
=
\exp\!\left(t_1x-\frac{t_2x^2}{2}+\frac{t_3x^3}{3}-\cdots\right)
=
\det(I_t+xM),
\]
\[
\sum_{n=0}^\infty p_n\frac{x^n}{n!}
=
\exp\!\left(t_1x+\frac{t_2x^2}{2}+\frac{t_3x^3}{3}+\cdots\right)
=
\det(I_t-xM)^{-1}.
\]
The corresponding recursions are
\[
a_{n+1}=\sum_{j=0}^n \binom{n}{j}(-1)^j j!\,a_{n-j}\,t_{j+1},
\qquad
p_{n+1}=\sum_{j=0}^n \binom{n}{j}j!\,p_{n-j}\,t_{j+1},
\]
with \(a_0=p_0=1\). Since \(\operatorname{rank}(G_n)\le t\), \(\det(G_n)=0\) for \(n>t\), hence \(a_n=0\) in that regime. The expected coefficients of the characteristic and permanental polynomials factor through these sequences via
\[
\mathbb E[b_i]=\binom{n}{i}a_i,\qquad
\mathbb E[d_i]=\binom{n}{i}p_i.
\]
This framework shows that, for independent sampled columns, average Gramian information can be compressed to the second-moment matrix \(M\) and its power traces, without assuming Gaussianity or independence of coordinates within a column [1309.2599].

## 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 \(A\in\mathbb R^{m\times n}\) with columns \(a_i\), one samples indices \(i_1,\dots,i_c\) independently with replacement according to probabilities \(p=\{p_i\}_{i=1}^n\) and forms
\[
\widehat G
=
\frac{1}{c}\sum_{t=1}^c \frac{1}{p_{i_t}}\, a_{i_t}a_{i_t}^T
=
(AS)(AS)^T.
\]
The estimator is unbiased:
\[
\mathbb E[\widehat G]=AA^T.
\]
This gives a literal average-over-samples approximation to the Gram matrix, with weights chosen to compensate for nonuniform sampling [1310.1502].

The same paper gives exact and probabilistic characterizations. If \(A=U\Sigma V^T\) has rank \(k\), selected columns \(a_{t_1},\dots,a_{t_c}\), and diagonal weights \(W=\operatorname{diag}(w_1,\dots,w_c)\), then
\[
AA^T=\sum_{j=1}^c w_j\,a_{t_j}a_{t_j}^T
\quad\text{if and only if}\quad
V_1^T W V_1 = I_k,
\]
equivalently \(V_1^T W^{1/2}\) has orthonormal rows. In the special case \(c=k=\operatorname{rank}(A)\), the unique diagonal weights are
\[
w_j=\frac{1}{\|V^T e_{t_j}\|_2^2}.
\]
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
\[
\operatorname{sr}(A)=\frac{\|A\|_F^2}{\|A\|_2^2}
\]
or on \(\operatorname{rank}(A)\), not on ambient matrix dimensions. With nearly-optimal probabilities \(p_i^\beta\ge \beta p_i^{\mathrm{opt}}\), \(0<\beta\le 1\), one obtains
\[
\frac{\|\widehat G-AA^T\|_2}{\|AA^T\|_2}\le \epsilon
\]
with probability at least \(1-\delta\), provided
\[
c \ge c_0(\epsilon)\,\frac{\operatorname{sr}(A)}{\beta\epsilon^2}
\ln\!\Big(\frac{\operatorname{rank}(A)}{\delta}\Big),
\qquad
c_0(\epsilon)=2+\frac{2\epsilon}{3}.
\]
An intrinsic-dimension variant replaces \(\ln(\operatorname{rank}(A)/\delta)\) by \(\ln(4\,\operatorname{sr}(A)/\delta)\), and leverage-score sampling yields a rank-based bound
\[
c \ge c_0(\epsilon)\,\frac{\operatorname{rank}(A)}{\epsilon^2}
\ln\!\Big(\frac{\operatorname{rank}(A)}{\delta}\Big).
\]
The practical guidance given is to prefer norm-squared probabilities
\[
p_i^{\mathrm{opt}}=\frac{\|a_i\|_2^2}{\|A\|_F^2},
\]
which minimize \(\mathbb E\|\widehat G-AA^T\|_F^2\) and empirically yield smaller spectral errors than leverage scores when \(\operatorname{sr}(A)\) is small [1310.1502].

## 7. Gramian averages in spherical designs and moment matching

For real and complex spherical designs, the Gramian is the matrix of pairwise inner products,
\[
\Gamma=[g_{ij}],\qquad g_{ij}=\langle x_i,x_j\rangle,
\]
and with weights \(w_i\) one may pass to the weighted Gram
\[
\Gamma_w = W^{1/2}\Gamma W^{1/2},\qquad W=\operatorname{diag}(w_1,\dots,w_N).
\]
If a weighted spherical design is defined as a cubature rule for a unitarily invariant polynomial space \(P\), then any unitary image is also such a design, and spherical designs for \(P\) are determined up to unitary equivalence by their Gramian. The reproducing kernel has the form
\[
K_P(x,y)=F(\langle x,y\rangle),
\]
and the associated potential
\[
A_{P,w}(X)=\sum_{i,j} w_i w_j\,F(\langle x_i,x_j\rangle)\ge 0
\]
vanishes exactly for weighted spherical designs. In this setting, Gramian-based energies encode averaged polynomial moments of pairwise inner products [2511.07452].

The paper makes the “average Gramian” viewpoint explicit through moment identities. For balanced sets,
\[
\sum_i w_i x_i = 0,
\]
and for the second moment,
\[
\sum_i w_i x_i x_i^\top = \frac{1}{d}I_d \quad \text{(real)},
\qquad
\sum_i w_i x_i x_i^* = \frac{1}{d}I_d \quad \text{(complex)}.
\]
Equivalently, the average outer product is isotropic. For unit weights \(w_i=1/N\) and \(X=[x_1\ \cdots\ x_N]\), this is the tight-frame identity
\[
XX^*=\frac{N}{d}I_d.
\]
In Gramian terms, \(\Gamma=X^*X\) has rank \(d\) with \(d\) nonzero eigenvalues all equal to \(N/d\), and the frame potential
\[
\operatorname{FP}(X)=\sum_{i,j}|\langle x_i,x_j\rangle|^2=\|\Gamma\|_F^2=\operatorname{Tr}(\Gamma^2)
\]
satisfies
\[
\operatorname{FP}(X)\ge \frac{N^2}{d},
\]
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
\[
F(t)=t^m-b_m(d)
\]
for real half-designs,
\[
F(z)=z^p\overline z^{\,q}-b_{p,q}(d)
\]
for complex \((p,q)\)-designs, and
\[
F(z)=|z|^{2t}-b_{t,t}(d)
\]
for complex projective \(t\)-designs. The design condition is that the sample Gramian reproduces spherical averages such as
\[
\sum_{i,j} w_i w_j \langle x_i,x_j\rangle^m = b_m(d).
\]
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 [2511.07452].

Source: https://www.emergentmind.com/topics/average-gramian-matrix