Gaussian Priors over Matrices and the Symmetric Kronecker Product

Summary

To build a concrete probabilistic linear solver we need tractable probability distributions over matrices. This note collects the required machinery: vectorisation and the Kronecker product (Ch. 15), the matrix-variate Gaussian, the Kronecker-structured covariance ${\Sigma }_0=V_0\otimes W_0$ that reduces inference cost from $\mathcal{O}(M^3)$-plus to $\mathcal{O}(N{M}^2)$, and the symmetric Kronecker product $W{\otimes }_{\ominus }W$ used to encode matrix symmetry as a linear constraint. These objects are the building blocks of the update equations in Probabilistic Linear Solvers - Algorithmic Scaffold.

Overview

A matrix $A\in {\mathbb{R}}^{N\times M}$ is treated here not as an algebraic operator but as a table of $NM$ numbers. The Gaussian family over vectors is transplanted onto matrices via vectorisation. The challenge is that a fully general covariance over $NM$ elements is $NM\times NM$ — larger than the matrix we are trying to invert — so structural restrictions are mandatory. Kronecker structure supplies the needed tractability without giving up support over all matrices.

Main Content

Vectorisation and the Kronecker product

Vectorisation ^def-vec

For $A\in {\mathbb{R}}^{N\times M}$, $\vec{A}\in {\mathbb{R}}^{NM}$ stacks the elements row after row (“row-major”), indexed by the same double index $(ij)$ as the matrix. The inverse operation is written $A=♮\vec{A}$.

Kronecker product ^def-kron

For $A\in {\mathbb{R}}^{N_A\times M_A}$, $B\in {\mathbb{R}}^{N_B\times M_B}$, the Kronecker product $A\otimes B$ has size $N_A{N}_B\times M_A{M}_B$ with elements $[A\otimes B{]}_{ij,k\ell }=[A{]}_{ik}[B{]}_{j\ell }$ (15.2). It “translates” between matrix multiplication and vectorisation:

$$[(A\otimes B)\vec{C}{]}_{ij}={\sum }_{k\ell }[A{]}_{ik}[C{]}_{k\ell }[B{]}_{j\ell }=[\overline{AC{B}^{\top }}{]}_{ij}.\text{\hspace{1em}(15.3)}$$

Key identities (all assuming the inverses/decompositions exist):

$$(A\otimes B)(C\otimes D)=AC\otimes BD,\text{\hspace{1em}(15.4)}$$ $$(A\otimes B{)}^{-1}=A^{-1}\otimes B^{-1},(A\otimes B{)}^{\top }=A^{\top }\otimes B^{\top }.\text{\hspace{1em}(15.5, 15.7)}$$

The matrix-variate Gaussian

Gaussian over matrix elements ^def-matrix-gaussian

For $X\in {\mathbb{R}}^{N\times K}$, define

$$\mathcal{N}(X;X_0,{\Sigma }_0):=\frac{\exp (-\frac{1}{2}(\vec{X}-\overrightarrow{X_0}{)}^{\top }{\Sigma }_0^{-1}(\vec{X}-\overrightarrow{X_0}))}{(2\pi {)}^{NK/2}|{\Sigma }_0{|}^{1/2}},\text{\hspace{1em}(19.3)}$$

where $\overrightarrow{X_0}\in {\mathbb{R}}^{NK}$ is a vectorised mean matrix and ${\Sigma }_0\in {\mathbb{R}}^{NK\times NK}$ is a symmetric positive-(semi)definite covariance. If ${\Sigma }_0$ is full rank (spd), this assigns non-zero density to every matrix — symmetric and asymmetric, definite and indefinite, invertible and non-invertible. (The arrow over the matrix argument is suppressed in shorthand.)

The Wishart distribution $\mathcal{W}(X;V,\nu )\propto |X{|}^{(\nu -N-1)/2}{e}^{-\frac{1}{2}\mathrm{tr}(V^{-1}X)}$ might seem natural (it lives on spd matrices, $X={\sum }_i{w}_i{w}_i^{\top }$ with $w_i\sim \mathcal{N}(0,V)$), but the posterior from conditioning on linear projections $Y=AS$ is not Wishart and has no compact form. Gaussians are used instead because linear projections keep them Gaussian.

Kronecker covariance for efficiency

A diagonal ${\Sigma }_0$ (independent elements) makes the Gram matrix cost $\mathcal{O}(N{M}^2)$ but its inverse costs more than $M^3$ — still too much. The fix is a Kronecker product covariance:

Kronecker covariance ^def-kron-cov

$${\Sigma }_0=V_0\otimes W_0,V_0,W_0\in {\mathbb{R}}^{N\times N}\text{spd}.\text{\hspace{1em}(19.8)}$$

Since a Kronecker product of two spd matrices is spd, this is a valid spd covariance assigning non-vanishing density to every matrix, while offering only $2\cdot \frac{1}{2}N(N+1)=N(N+1)$ degrees of freedom (versus $\frac{1}{2}{N}^2(N^2+1)$ for a general ${\Sigma }_0$). Element covariances factorise:

$$\mathrm{cov}(A_{ij},A_{k\ell })=[V_0\otimes W_0{]}_{ij,k\ell }=[V_0{]}_{ik}[W_0{]}_{j\ell },\text{\hspace{1em}(19.9)}$$

and marginal variances depend only on the diagonals: $\mathrm{var}(A_{ij})=[V_0{]}_{ii}[W_0{]}_{jj}$.

Caution (Fig. 19.2): it is tempting but not entirely correct to read $V_0$ as “covariance among rows” and $W_0$ “among columns”. A matrix $\tilde{A}=B\odot C+A_0$ built from column-vectors $b_i\sim \mathcal{N}(0,V_0)$ and row-vectors $c_j\sim \mathcal{N}(0,W_0)$ satisfies (19.9) but is not Gaussian. The Kronecker restriction limits complexity, not the support — the measure still covers all matrices.

The symmetric Kronecker product (encoding symmetry)

Symmetry $A_{ij}-A_{ji}=0$ is a linear constraint, and Gaussians are closed under linear constraints (positive-definiteness, being a cone, is not — see ^cone). Two orthogonal projection operators on ${\mathbb{R}}^{N^2}$ (space of vectorised $N\times N$ matrices) formalise this:

Symmetric / anti-symmetric projections ^def-proj

$$[{\Pi }_{\ominus }{]}_{(ij),(k\ell )}:=\frac{1}{2}({\delta }_{ik}{\delta }_{j\ell }+{\delta }_{i\ell }{\delta }_{jk}),♮{\Pi }_{\ominus }\vec{X}=\frac{1}{2}(X+X^{\top }),$$ $$[{\Pi }_{\oplus }{]}_{(ij),(k\ell )}:=\frac{1}{2}({\delta }_{ik}{\delta }_{j\ell }-{\delta }_{i\ell }{\delta }_{jk}),♮{\Pi }_{\oplus }\vec{X}=\frac{1}{2}(X-X^{\top }).$$

They are orthogonal projectors that jointly span ${\mathbb{R}}^{N^2}$: ${\Pi }_{\ominus }{\Pi }_{\ominus }={\Pi }_{\ominus }$, ${\Pi }_{\oplus }{\Pi }_{\oplus }={\Pi }_{\oplus }$, ${\Pi }_{\ominus }{\Pi }_{\oplus }=0$, ${\Pi }_{\ominus }+{\Pi }_{\oplus }=I_{N^2}$ (19.15).

A Kronecker product $W\otimes W$ splits into a symmetric and skew-symmetric part:

$$W\otimes W=\underset{=:W{\otimes }_{\ominus }W}{\underbrace{{\Pi }_{\ominus }(W\otimes W){\Pi }_{\ominus }}}+\underset{=:W{\otimes }_{\oplus }W}{\underbrace{{\Pi }_{\oplus }(W\otimes W){\Pi }_{\oplus }}}.$$

Symmetric / skew Kronecker product ^def-symkron

$$[C{\otimes }_{\ominus }D{]}_{(ij),(k\ell )}=\frac{1}{2}(C_{ik}{D}_{j\ell }+C_{i\ell }{D}_{jk}),[C{\otimes }_{\oplus }D{]}_{(ij),(k\ell )}=\frac{1}{2}(C_{ik}{D}_{j\ell }-C_{i\ell }{D}_{jk}).\text{\hspace{1em}(19.16)}$$

They inherit some Kronecker properties. In particular, on symmetric matrices $(W{\otimes }_{\ominus }W{)}^{-1}=W^{-1}{\otimes }_{\ominus }{W}^{-1}$, but in general $(C{\otimes }_{\oplus }D{)}^{-1}\ne C^{-1}{\otimes }_{\oplus }{D}^{-1}$. Also $(C{\otimes }_{\ominus }D)\vec{X}=\frac{1}{2}(CX{D}^{\top }+C{X}^{\top }{D}^{\top })$ (19.17).

Symmetry information is written as a Dirac likelihood $p(\ominus \mid A)=\delta ({\Pi }_{\oplus }\vec{A}-0)={\lim }_{\beta \to 0}\mathcal{N}(\vec{0};{\Pi }_{\oplus }\vec{A},\beta I)$ (19.18). Conditioning the Gaussian prior on it yields a posterior that, with Kronecker covariance ${\Sigma }_0=W_0\otimes W_0$, collapses to a symmetric prior mean and a symmetric Kronecker covariance:

$$p(A)=\mathcal{N}(A;A_0,W_0{\otimes }_{\ominus }{W}_0),\text{\hspace{1em}(19.20)}$$

which assigns non-zero measure only to symmetric matrices.

Examples

Degrees of freedom

A full spd covariance on $N\times N$ matrix elements has $\frac{1}{2}{N}^2(N^2+1)$ parameters. Kronecker $V_0\otimes W_0$ has $N(N+1)$. The symmetric $N\times N$ matrices form a space of dimension $\frac{1}{2}N(N+1)$; if $W$ is full rank, $W{\otimes }_{\ominus }W$ has rank $\frac{1}{2}N(N+1)$ and its inverse on that space is $W^{-1}{\otimes }_{\ominus }{W}^{-1}$.

Sampling intuition (Figs. 19.3–19.4)

Samples from $\mathcal{N}(I,V_0\otimes W_0)$ with $V_0=W_0=I$ look like uncorrelated noise; with $V_0=\mathrm{diag}[10^4,9^4,8^4,\dots ]$ the row structure dominates. Samples from the symmetry-encoding $\mathcal{N}(A_0,W_0{\otimes }_{\ominus }{W}_0)$ are visibly symmetric matrices.

Connections

• Provides the covariance objects used in the mean/covariance update equations of Probabilistic Linear Solvers - Algorithmic Scaffold.
• The matrix inversion lemma and block-inverse identities (Ch. 15.5) are the algebraic tools for turning a posterior over $A$ into an estimator $H_i=A_M^{-1}$.
• Builds directly on Gaussian Distributions and Algebra (closure under linear maps and conditioning).

See Also

• Gaussian Distributions and Algebra — the vector-Gaussian identities transplanted here.
• Probabilistic Linear Solvers - Algorithmic Scaffold — uses these priors to instantiate the solver.
• The Linear Algebra Problem and Evaluation Strategies — the observation model $Y=AS$ conditioned on.