Can the covariance matrix have imaginary eigenvalues?

The covariance matrix

The covariance matrix is commonly defined as:

\[ \begin{aligned} &\text{cov}(X,Y) := \mathbb{E}[(X - \mathbb{E}(X))(Y - \mathbb{E}(Y))^T] \\ &\text{where} \\ &\quad\quad X, Y : (\Omega, \mathcal{F}, P) \xrightarrow{\text{measurable}} (\mathbb{R}^n, \mathcal{B}(\mathbb{R}^n)) \\ &\quad\quad (\Omega, \mathcal{F}, P) \text{ is a probability space} \\ &\quad\quad \mathcal{B}(\mathbb{R}^n) \text{ is the Borel sigma-algebra on } \mathbb{R}^n \end{aligned} \]

If we are given only one random variable \(X\), we can define the covariance matrix of \(X\) as \(\text{cov}(X) := \text{cov}(X,X)\), which is also called the auto-covariance matrix of \(X\). Clearly this auto-covariance matrix is real and symmetric:

\[ \begin{aligned} \text{cov}(X) &= \mathbb{E}[(X - \mathbb{E}(X))(X - \mathbb{E}(X))^T] \\ &= \begin{bmatrix} \mathbb{E}[(X_1 - \mathbb{E}(X_1))^2] & \cdots & \mathbb{E}[(X_1 - \mathbb{E}(X_1))(X_n - \mathbb{E}(X_n))] \\ \vdots & \ddots & \vdots \\ \mathbb{E}[(X_n - \mathbb{E}(X_n))(X_1 - \mathbb{E}(X_1))] & \cdots & \mathbb{E}[(X_n - \mathbb{E}(X_n))^2] \end{bmatrix} \\ &= \text{cov}(X)^T \end{aligned} \]

As a real symmetric matrix, it cannot have imaginary eigenvalues. Remember that an eigenvalue \(\lambda\) of a matrix \(A\) is defined by the equation

\[ Av = \lambda v, \quad v \neq 0, \]

where we allow \(\lambda \in \mathbb{C}\) and \(v \in \mathbb{C}^n\), since \(\mathbb{C}\) is where the characteristic polynomial \(\chi_A(\lambda) = \det(A - \lambda I)\) splits. The claim is that for symmetric real \(A\) we never leave \(\mathbb{R}\).

NoteWhy the eigenvalues of a real symmetric matrix are real

Let \(Av = \lambda v\) with \(v \in \mathbb{C}^n \setminus \{0\}\). Then

\[ \lambda \, \bar{v}^T v = \bar{v}^T A v \overset{(1)}{=} (A \bar{v})^T v \overset{(2)}{=} \overline{(A v)}^T v = \overline{(\lambda v)}^T v = \bar{\lambda} \, \bar{v}^T v, \]

where \((1)\) uses \(A^T = A\) and \((2)\) uses \(\bar{A} = A\). Since \(\bar{v}^T v = \sum_i |v_i|^2 > 0\) we may divide by it and obtain \(\lambda = \bar{\lambda}\), i.e. \(\lambda \in \mathbb{R}\). \(\square\)

It is worth noting exactly which two hypotheses that argument spends: \(A^T = A\) and \(\bar A = A\). Together they say

\[ A^* = A, \quad \text{where } A^* := \bar{A}^T, \]

i.e. \(A\) is fixed by the involution \((-)^*\) on matrices. Reality of the spectrum is a statement about the involution, not about symmetry alone. Keep an eye on this — it is the hypothesis we are about to break, even though we will set out to break a different one.

Breaking the symmetry

How would we be able to have a non-symmetric auto-covariance matrix? We would need to have

\[ \exists i, j \in \{1, \ldots, n\} : \mathbb{E}[(X_i - \mathbb{E}(X_i))(X_j - \mathbb{E}(X_j))] \neq \mathbb{E}[(X_j - \mathbb{E}(X_j))(X_i - \mathbb{E}(X_i))]. \]

Here, let us first establish without loss of generality that all \(X_i\) are centered, i.e. \(\mathbb{E}(X_i) = 0\) for all \(i\). Then the above condition can be rewritten as:

\[ \exists i, j \in \{1, \ldots, n\} : \mathbb{E}[X_i X_j] \neq \mathbb{E}[X_j X_i]. \]

The multiplication of real random variables is commutative, as it is lifted from the multiplication of real numbers. We will have to consider a non-commutative multiplication to break the symmetry of the auto-covariance matrix, obviously! We can now move over to another ring. In this example, we will consider \(\mathbb{H}\), the ring of quaternions, which is a non-commutative extension of the complex numbers. The quaternions are defined as:

\[ \mathbb{H} := \{a + bi + cj + dk : a, b, c, d \in \mathbb{R}, i^2 = j^2 = k^2 = ijk = -1\} \]

Alternatively, we could mod out the following ideal from the polynomial ring \(\mathbb{R}[i,j,k]\):

\[ \mathbb{H} \cong_\text{Ring} \mathbb{R}[i,j,k] / (i^2 + 1, j^2 + 1, k^2 + 1, ij - k, jk - i, ki - j). \]

We call the elements \(bi + cj + dk\) imaginary, in analogy to the complex numbers. The elements \(a\) are called real. Conjugation reverses the sign of the imaginary part, \(\overline{a + bi + cj + dk} = a - bi - cj - dk\), and — this will matter — it is an anti-automorphism: \(\overline{pq} = \bar q \bar p\).

A first attempt: take the constant random variable

\[ X = \begin{bmatrix} i \\ j \end{bmatrix} \in \mathbb{H}^2 \]

and use \(ij = k\), \(ji = -k\). Its auto-covariance matrix is

\[ \text{cov}(X) = \begin{bmatrix} ii & ij \\ ji & jj \end{bmatrix} = \begin{bmatrix} -1 & k \\ -k & -1 \end{bmatrix} \neq \text{cov}(X)^T. \]

Non-symmetric — the plan worked. So what are its eigenvalues?

Eigenvalues over a non-commutative ring

Here everything we know breaks at once, and it is worth being careful, because this is where I originally fooled myself.

There are two kinds of eigenvalue. Over \(\mathbb{H}\) the equations \(Av = \lambda v\) and \(Av = v\lambda\) are genuinely different, and they define the left and right eigenvalues of \(A\). Left eigenvalues are badly behaved: a \(2 \times 2\) quaternionic matrix can have infinitely many of them, and there is no useful spectral theory (Wood 1985). Right eigenvalues are the right notion, and from here on “eigenvalue” means right eigenvalue.

Eigenvalues come in spheres. If \(Av = v\lambda\) and \(u \in \mathbb{H}^\times\), then \(A(vu) = (vu)(u^{-1}\lambda u)\), so the whole similarity class \([\lambda] = \{u^{-1}\lambda u : u \in \mathbb{H}^\times\}\) consists of eigenvalues. For a non-real \(\lambda = a + b\,\mathrm{i}\) this class is the 2-sphere \(\{a + b(xi + yj + zk) : x^2+y^2+z^2 = 1\}\) of radius \(|b|\) around \(a\). This is the honest version of a remark that is easy to get wrong: the solutions of \(\lambda^2 = -1\) in \(\mathbb{H}\) are not \(\{\pm i, \pm j, \pm k\}\), they are all of \(S^2 \subset \operatorname{Im}\mathbb{H}\). Each similarity class has exactly one representative in \(\mathbb{C}\) with non-negative imaginary part, and that is the number we will quote.

There is no determinant. A determinant needs commutativity to be multiplicative; over \(\mathbb{H}\) the best available substitute is the Dieudonné determinant, which takes values in \(\mathbb{R}_{\geq 0}\) and cannot detect a phase. So \(\det(A - \lambda I) = 0\) is not a valid way to find eigenvalues, and any computation that uses it is a computation with an undefined symbol.

What works instead. Every quaternion splits uniquely as \(q = z_1 + z_2 j\) with \(z_1, z_2 \in \mathbb{C}\), and multiplication is faithfully recorded by the complex adjoint: writing \(A = A_1 + A_2 j\) with \(A_1, A_2 \in \mathbb{C}^{n \times n}\), set

\[ \chi(A) := \begin{bmatrix} A_1 & A_2 \\ -\overline{A_2} & \overline{A_1} \end{bmatrix} \in \mathbb{C}^{2n \times 2n}. \]

Then \(\chi\) is an injective homomorphism of real algebras, the \(2n\) eigenvalues of \(\chi(A)\) come in conjugate pairs, and the \(n\) of them with non-negative imaginary part are exactly the standard representatives of the right eigenvalues of \(A\) (Zhang 1997). In other words: quaternionic linear algebra is complex linear algebra with a symmetry constraint, and numpy can do all of it.

import numpy as np

def adjoint(Q):
    """Complex adjoint of a quaternionic matrix.

    `Q` has shape (n, n, 4): the coefficients in the basis (1, i, j, k).
    """
    A1 = Q[..., 0] + 1j * Q[..., 1]   # the 1 and i part
    A2 = Q[..., 2] + 1j * Q[..., 3]   # the j and k part, since c j + d k = (c + d i) j
    return np.block([[A1, A2], [-A2.conj(), A1.conj()]])

def right_eigenvalues(Q):
    """Standard representatives of the right eigenvalues of a quaternionic matrix."""
    ev = np.linalg.eigvals(adjoint(Q))
    ev = ev[np.lexsort((ev.imag, ev.real))]   # sort so conjugate pairs are adjacent
    return ev[1::2]                           # keep the one with imaginary part >= 0

# cov(X) for X = [i, j]^T
C = np.array([[(-1, 0, 0, 0), (0, 0, 0, 1)],
              [(0, 0, 0, -1), (-1, 0, 0, 0)]], dtype=float)

print(np.round(right_eigenvalues(C), 6))

[-2.+0.j -0.+0.j]

Real. Both of them.

And this is not bad luck, it is a theorem. Every component of \(X = [i, j]^T\) is a purely imaginary quaternion, so \(\bar X_m = -X_m\), and because conjugation reverses products,

\[ \overline{X_m X_n} = \bar X_n \bar X_m = (-X_n)(-X_m) = X_n X_m . \]

Taking expectations, \(\overline{\text{cov}(X)_{mn}} = \text{cov}(X)_{nm}\), which is to say

\[ \text{cov}(X)^* = \text{cov}(X). \]

The matrix is not symmetric, but it is Hermitian, and Hermitian quaternionic matrices have real eigenvalues for exactly the reason real symmetric matrices do. We broke commutativity and the involution quietly patched the hole.

A covariance matrix that really does have imaginary eigenvalues

To get past this we have to stop the conjugation from cancelling, which means the components of \(X\) must mix real and imaginary parts. Let \(\varepsilon\) take the values \(+1\) and \(-1\) with probability \(\tfrac12\) each, and set

\[ X := \varepsilon \begin{bmatrix} 1 + i \\ j \end{bmatrix} \in \mathbb{H}^2 . \]

This one is an honest random variable and it is centered, \(\mathbb{E}(X) = 0\), so no generality was lost above. Since \(\varepsilon^2 = 1\),

\[ \text{cov}(X) = \begin{bmatrix} (1+i)(1+i) & (1+i)j \\ j(1+i) & jj \end{bmatrix} = \begin{bmatrix} 2i & j + k \\ j - k & -1 \end{bmatrix}, \]

which is neither symmetric nor Hermitian.

C = np.array([[(0, 2, 0, 0), (0, 0, 1, 1)],
              [(0, 0, 1, -1), (-1, 0, 0, 0)]], dtype=float)

print(np.round(right_eigenvalues(C), 6))

[-1.+2.j  0.+0.j]

There it is: the eigenvalues are \(0\) and \(-1 + 2i\), and by the discussion above the second one is really the entire sphere of radius \(2\) centered at \(-1\) in \(\mathbb{H}\). So yes — a covariance matrix can have imaginary eigenvalues.

The uncomfortable part is what it cost us, and we will come back to that at the end. First, what would imaginary eigenvalues even mean for data?

Rotation of data

To find out the geometric significance of the constructed imaginary eigenvalues, we first need to understand in what way the eigenvalues and eigenvectors of the auto-covariance matrix encode geometric properties of the data.

Principal components as eigenvectors

When working with collected samples \(x_1, \ldots, x_N\) of a random variable \(X\), we often work with the sample covariance matrix \(\Sigma_X\), which is defined using the data matrix \(D_X\):

\[ \begin{aligned} &D_X := \begin{bmatrix} x_1^T \\ x_2^T \\ \vdots \\ x_N^T \end{bmatrix} \in \mathbb{R}^{N \times n} \end{aligned} \]

The sample covariance matrix is then defined as:

\[ \begin{aligned} &\Sigma_X := \frac{1}{N} D_X^T D_X \\ &= \frac{1}{N} \sum_{i=1}^N x_i x_i^T \end{aligned} \]

It is well known (see corollary 5.2 https://www.cs.columbia.edu/~djhsu/AML/lectures/notes-pca.pdf), that the eigenvectors of \(\Sigma_X\) are the principal components of the data, and that the eigenvalues of \(\Sigma_X\) are the variances of the data along those principal components. The principal components are obtained by rotating the data in such a way that the variance along each principal component is maximized:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

rng = np.random.default_rng(42)

angle = np.radians(40)
R = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
X = rng.multivariate_normal([1, 3], R @ np.diag([9, 1]) @ R.T, size=1000)

eigenvalues, eigenvectors = np.linalg.eigh(np.cov(X.T))
eigenvalues, eigenvectors = eigenvalues[::-1], eigenvectors[:, ::-1]
mean = X.mean(axis=0)
colors = ["#e63946", "#457b9d"]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.scatter(X[:, 0], X[:, 1], alpha=0.25, s=8, color="gray")
for i, (lam, vec, c) in enumerate(zip(eigenvalues, eigenvectors.T, colors)):
    std = np.sqrt(lam)
    arrow_kw = dict(arrowstyle="-|>", lw=2, color=c)
    ax1.annotate("", xy=mean + std * vec, xytext=mean, arrowprops=arrow_kw)
    ax1.annotate("", xy=mean - std * vec, xytext=mean, arrowprops=arrow_kw)
    ax1.text(*(mean + (std + 0.4) * vec), rf"$\lambda_{i+1}={lam:.1f}$, $\sigma={std:.1f}$",
             color=c, fontsize=9, ha="center")
ax1.set_aspect("equal")
ax1.grid(True, linestyle=":", alpha=0.5)
ax1.set_title("Eigenvectors as principal components\n(arrow length = std along PC)")

projections = (X - mean) @ eigenvectors
for i, (lam, c) in enumerate(zip(eigenvalues, colors)):
    proj = projections[:, i]
    xs = np.linspace(proj.min(), proj.max(), 300)
    ax2.hist(proj, bins=50, alpha=0.55, color=c, density=True,
             label=rf"PC {i+1}: var={proj.var():.2f}, $\lambda$={lam:.2f}")
    ax2.plot(xs, norm.pdf(xs, scale=np.sqrt(lam)), color=c, lw=2)
ax2.set(xlabel="Projection value", ylabel="Density",
        title="Variance along each PC = eigenvalue")
ax2.legend(fontsize=9)
ax2.grid(True, linestyle=":", alpha=0.5)

plt.tight_layout()
plt.show()

Here, the data matrix would look like this:

\[ D_X = \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \end{bmatrix} \]

and the sample covariance matrix would be:

\[ \Sigma_X = \frac{1}{N} D_X^T D_X = \frac{1}{N} \begin{bmatrix} x_{11} & x_{21} & \cdots & x_{N1} \\ x_{12} & x_{22} & \cdots & x_{N2} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{N1} & x_{N2} \end{bmatrix} = \frac{1}{N} \begin{bmatrix} \sum_{i=1}^N x_{i1}^2 & \sum_{i=1}^N x_{i1} x_{i2} \\ \sum_{i=1}^N x_{i1} x_{i2} & \sum_{i=1}^N x_{i2}^2 \end{bmatrix} \]

Note

Remember that the standard inner product on \(\mathbb{R}^n\) is defined as \(\langle x, y \rangle = x^T y\), and we have the well known cosine similarity formula: \[ \cos(\theta)\|x\| \|y\| = \langle x, y \rangle = x^T y = \sum_{i=1}^n x_i y_i \] If \(x_{\cdot ,1}\) and \(x_{\cdot, 2}\) in the above example were normed, then the covariance matrix would encode the cosine similarity between the rows of the data matrix, or equivalently the sampled \(x_i\).

The diagonals would clearly just be the lengths (which we assumed to be normed). Dropping our norm assumption, we would therefore be able to obtain the cosine similarity of the data points, by considering the off-diagonal entries of the covariance matrix and dividing them by the square root of the product of the diagonal entries.

Similarly, we can detect the principal components of 3d point clouds. This however crucially only detects the spread of the data along the principal components. If our data has a shape different from a Gaussian blob, PCA will fall short.

import numpy as np
import plotly.graph_objects as go

rng = np.random.default_rng(42)
COLORS = ["#e63946", "#457b9d", "#2a9d8f"]

def add_pca_arrows(fig, X):
    vals, vecs = np.linalg.eigh(np.cov(X.T))
    vals, vecs = vals[::-1], vecs[:, ::-1]
    mean = X.mean(axis=0)
    for i, (lam, vec, c) in enumerate(zip(vals, vecs.T, COLORS)):
        scale = np.sqrt(lam)
        for sign in [1, -1]:
            end = mean + sign * scale * vec
            fig.add_trace(go.Scatter3d(
                x=[mean[0], end[0]], y=[mean[1], end[1]], z=[mean[2], end[2]],
                mode="lines", line=dict(color=c, width=6),
                showlegend=(sign == 1), name=f"PC {i+1} (λ={lam:.2f})"
            ))
            fig.add_trace(go.Cone(
                x=[end[0]], y=[end[1]], z=[end[2]],
                u=[sign * 0.3 * vec[0]], v=[sign * 0.3 * vec[1]], w=[sign * 0.3 * vec[2]],
                colorscale=[[0, c], [1, c]], showscale=False, showlegend=False,
                sizemode="absolute", sizeref=0.4
            ))

layout_kw = dict(margin=dict(l=0, r=0, t=40, b=0), scene=dict(aspectmode="data"))

# Plot 1: 3D Gaussian blob
cov3 = np.array([[4, 2, 0.5], [2, 3, 1], [0.5, 1, 1]])
X1 = rng.multivariate_normal([0, 0, 0], cov3, size=800)

fig1 = go.Figure(layout=layout_kw)
fig1.add_trace(go.Scatter3d(x=X1[:,0], y=X1[:,1], z=X1[:,2], mode="markers",
    marker=dict(size=2, color="#555", opacity=0.5), name="samples"))
add_pca_arrows(fig1, X1)
fig1.update_layout(title="3D Gaussian blob with principal components")
fig1.show()

# Plot 2: xz Gaussian blob with bell-curve bump in y
xz = rng.multivariate_normal([0, 0], [[4, 1], [1, 2]], size=800)
x2, z2 = xz[:, 0], xz[:, 1]
y2 = 2.5 * np.exp(-(x2**2 + z2**2) / 5) + rng.normal(0, 0.1, len(x2))
X2 = np.column_stack([x2, y2, z2])

fig2 = go.Figure(layout=layout_kw)
fig2.add_trace(go.Scatter3d(x=X2[:,0], y=X2[:,1], z=X2[:,2], mode="markers",
    marker=dict(size=2, color="#555", opacity=0.5), name="samples"))
add_pca_arrows(fig2, X2)
fig2.update_layout(title="xz Gaussian blob with bell-curve bump along y: falls short")
fig2.show()

Some standard solutions to this problem are Nonlinear PCA, for example by using kernel methods.

None of these approaches take into account the orientation of our data. It does not matter in which order the data points \(x_i\) are collected and assembled into the data matrix \(D_X\): Assume we permute the data points, resulting in a permuted data matrix \(P D_X\), where \(P\) is a permutation matrix. Then we have \[ \Sigma_{P X} = \frac{1}{N} (P D_X)^T (P D_X) = \frac{1}{N} D_X^T P^T P D_X = \frac{1}{N} D_X^T D_X = \Sigma_X. \]

What the imaginary part measures

Identify \(\mathbb{R}^3\) with the imaginary quaternions,

\[ \begin{aligned} \iota &: \mathbb{R}^3 \xrightarrow{\ \sim\ } \operatorname{Im}\mathbb{H}, \\ \iota(x,y,z) &= xi + yj + zk . \end{aligned} \]

Under this identification, quaternion multiplication is the dot product and the cross product at the same time: for \(u, v \in \operatorname{Im}\mathbb{H}\),

\[ u v = -\langle u, v \rangle + u \times v , \]

with the real part on the left and the imaginary part on the right. So if \(X = (X_1, \ldots, X_n)\) has purely imaginary components, its auto-covariance matrix splits as

\[ \text{cov}(X)_{mn} = \mathbb{E}[X_m X_n] = \underbrace{-\,\mathbb{E}\langle X_m, X_n \rangle}_{\text{real, symmetric}} \;+\; \underbrace{\mathbb{E}[X_m \times X_n]}_{\text{imaginary, antisymmetric}} \;=\; -\Sigma_{mn} + \Omega_{mn} . \]

The real part is just (minus) the ordinary covariance matrix. The imaginary part \(\Omega\) is the new object: the expected cross product, i.e. the mean oriented area — the angular momentum — swept out between the two variables. And the two halves are exactly what made the matrix Hermitian a moment ago: a real symmetric part plus an imaginary antisymmetric part is precisely what \(A^* = A\) says.

So the spectrum stays real. But that does not mean \(\Omega\) is empty. Recall \(\Sigma_{PX} = \Sigma_X\) from above: the ordinary covariance matrix cannot see the order of the samples. \(\Omega\) can, and it is the only part that can.

Take a noisy helix, sampled along its path, and compare it with the very same point cloud traversed backwards. PCA is identical — it is literally the same matrix — while \(\Omega\) flips sign and points along the axis of rotation:

t = np.linspace(0, 6 * np.pi, 600)
helix = np.stack([np.cos(t), np.sin(t), 0.25 * t], axis=1)
helix = helix + rng.normal(0, 0.03, helix.shape)
helix = helix - helix.mean(axis=0)
backwards = helix[::-1]

def Sigma(X):
    """The ordinary (sample) covariance matrix: the real, symmetric part."""
    return np.cov(X.T, bias=True)

def Omega(X):
    """The mean cross product of consecutive samples: the imaginary, antisymmetric part."""
    return np.cross(X[:-1], X[1:]).mean(axis=0)

assert np.allclose(Sigma(helix), Sigma(backwards))       # PCA cannot tell them apart
assert np.allclose(Omega(helix), -Omega(backwards))      # Omega can

fig3 = go.Figure(layout=layout_kw)
fig3.add_trace(go.Scatter3d(x=helix[:,0], y=helix[:,1], z=helix[:,2], mode="markers",
    marker=dict(size=2, color="#555", opacity=0.5), name="samples"))
add_pca_arrows(fig3, helix)

for X, label, c in [(helix, "as sampled", "#8338ec"), (backwards, "reversed", "#fb8500")]:
    w = Omega(X)
    end = 2.5 * w / np.linalg.norm(w)
    fig3.add_trace(go.Scatter3d(x=[0, end[0]], y=[0, end[1]], z=[0, end[2]],
        mode="lines", line=dict(color=c, width=8),
        name=f"Ω, {label} (|Ω|={np.linalg.norm(w):.3f})"))
    fig3.add_trace(go.Cone(x=[end[0]], y=[end[1]], z=[end[2]],
        u=[0.3*w[0]], v=[0.3*w[1]], w=[0.3*w[2]], colorscale=[[0, c], [1, c]],
        showscale=False, showlegend=False, sizemode="absolute", sizeref=0.4))

fig3.update_layout(title="Same cloud, same principal components, opposite Ω")
fig3.show()

The three principal components describe the box the helix lives in and nothing else. \(\Omega\) is a single vector along the helix axis whose sign is the handedness of the winding — the rotation that the friend at work had in mind, sitting in the imaginary part of the covariance matrix all along.

What we actually broke

None of this was really about quaternions. The general setting is a non-commutative probability space (Nica and Speicher 2006): a unital \(*\)-algebra \(\mathcal{A}\) together with a state \(\varphi : \mathcal{A} \to \mathbb{C}\), i.e. a linear functional with \(\varphi(1) = 1\) and \(\varphi(a^*a) \geq 0\). The elements of \(\mathcal{A}\) are the random variables and \(\varphi\) is the expectation; classical probability is the commutative case \(\mathcal{A} = L^\infty(\Omega, \mathcal{F}, P)\), \(\varphi = \mathbb{E}\).

In that language the covariance matrix of \(a_1, \ldots, a_n \in \mathcal{A}\) is \([\varphi(a_m^* a_n)]_{mn}\): the Gram matrix of the sesquilinear form \(\langle a, b \rangle_\varphi := \varphi(a^* b)\). Positivity of the state says this form is positive semi-definite, and that is the entire content — every Gram matrix is Hermitian with real, non-negative eigenvalues, and completing \(\mathcal{A}\) with respect to \(\langle -, - \rangle_\varphi\) is the GNS construction. No amount of non-commutativity will ever produce an imaginary variance. The example that worked escaped only because it was built from \(X X^T\) rather than \(X X^*\), and it paid for the escape with positivity: the matrix stopped being a Gram matrix, and its eigenvalues stopped being variances. What breaks is not commutativity. It is the involution.

What non-commutativity does break is \(\varphi(ab) = \varphi(ba)\). Splitting \(\varphi(a_m a_n)\) into symmetric and antisymmetric parts gives the ordinary covariance plus \(\tfrac{1}{2}\varphi([a_m, a_n])\), and the second term vanishes for all \(m, n\) exactly when \(\varphi\) is a trace. That is \(\Omega\) again, in general form. In quantum mechanics it is the canonical commutation relation \([q, p] = i\hbar\): the covariance matrix of a state is \(\sigma + \tfrac{i\hbar}{2}\Omega\) with \(\Omega\) the symplectic form, positive semi-definiteness of that combination is the Robertson–Schrödinger uncertainty relation, and Williamson’s theorem diagonalizes it not by eigenvalues but by symplectic eigenvalues — the generators of rotation in phase space (Adesso, Ragy, and Lee 2014). The rotation was real. It just lives in the imaginary part rather than in the spectrum.

Whether a trace exists at all is a structural property of \(\mathcal{A}\), and this is where von Neumann algebras answer the question we started with. A finite von Neumann algebra — a type \(\mathrm{II}_1\) factor such as the free group factor \(L(\mathbb{F}_n)\), the natural home of free probability (Voiculescu, Dykema, and Nica 1992) — carries a faithful normal tracial state, so its covariance matrices are symmetric again and we are back where we began. Type \(\mathrm{III}\) factors admit no trace whatsoever, and Tomita–Takesaki theory measures the failure by a canonical one-parameter group of automorphisms, the modular flow \(\sigma_t = \Delta^{it}(-)\Delta^{-it}\), which is trivial precisely when \(\varphi\) is tracial (Takesaki 2003); Connes’ classification of type \(\mathrm{III}\) factors is a classification by the spectrum of that flow (Connes 1973). So “can the covariance matrix be non-symmetric?” becomes “is the state a trace?”, and the rotation hiding in the imaginary part becomes, quite literally, a flow.

So: I was wrong, and I was also not. A covariance matrix in the sense of a Gram matrix of a state is Hermitian and has real, non-negative eigenvalues, no matter how non-commutative the underlying ring is — I was wrong for a better reason than the one I gave. But the rotation my friend and I were arguing about is genuinely there, in the antisymmetric imaginary part that classical statistics throws away before it ever writes the matrix down.

References

Adesso, Gerardo, Sammy Ragy, and Antony R. Lee. 2014. “Continuous Variable Quantum Information: Gaussian States and Beyond.” Open Systems & Information Dynamics 21 (01n02): 1440001. https://doi.org/10.1142/S1230161214400010.

Connes, Alain. 1973. “Une Classification Des Facteurs de Type III.” Annales Scientifiques de l’École Normale Supérieure 6 (2): 133–252.

Nica, Alexandru, and Roland Speicher. 2006. Lectures on the Combinatorics of Free Probability. Vol. 335. London Mathematical Society Lecture Note Series. Cambridge University Press.

Takesaki, Masamichi. 2003. Theory of Operator Algebras II. Vol. 125. Encyclopaedia of Mathematical Sciences. Springer.

Voiculescu, Dan V., Kenneth J. Dykema, and Alexandru Nica. 1992. Free Random Variables. Vol. 1. CRM Monograph Series. American Mathematical Society.

Wood, R. M. W. 1985. “Quaternionic Eigenvalues.” Bulletin of the London Mathematical Society 17 (2): 137–38.

Zhang, Fuzhen. 1997. “Quaternions and Matrices of Quaternions.” Linear Algebra and Its Applications 251: 21–57. https://doi.org/10.1016/0024-3795(95)00543-9.