Proving Cochran's Theorem from Orthogonal Subspaces

From Matrices to Subspaces

Cochran's theorem is often used to justify splitting a Gaussian sum of squares into components that are independent and have chi-square distributions. The underlying reason is geometric: each matrix selects an orthogonal subspace, and a complete collection of these subspaces supplies disjoint groups of Gaussian coordinates. The matrix ranks count the dimensions of those groups.

We develop the centered, spherical Gaussian case. The geometry comes first; once the projected subspaces are known to be mutually orthogonal and to fill the ambient space, the distributional result follows by choosing coordinates adapted to them.

Definition: A matrix \(P\in\mathbb{R}^{n\times n}\) is an orthogonal projection matrix if \(P=P'\) and \(P^2=P\). Its range is \(\mathcal{R}(P)=\{Px:x\in\mathbb{R}^n\}\), and its rank is the dimension of that range. For subspaces \(S_i,S_j\subseteq\mathbb{R}^n\), orthogonality means \(u'v=0\) for every \(u\in S_i\) and \(v\in S_j\).

Symmetry and idempotence are both essential to the projection interpretation. Idempotence says that applying \(P\) twice gives the same result as applying it once. Symmetry ensures that the part removed by the projection is perpendicular to the part retained.

Lemma: If \(P=P'\) and \(P^2=P\), then \(P\) is the orthogonal projection onto \(\mathcal{R}(P)\), and \(\ker(P)=\mathcal{R}(P)^\perp\). In particular, for every \(x\), \(x'Px=\|Px\|^2\).

To verify these claims, write \(x=Px+(I-P)x\). The first term belongs to \(\mathcal{R}(P)\), while \(P(I-P)x=(P-P^2)x=0\), so the second belongs to \(\ker(P)\). If \(u=Pa\in\mathcal{R}(P)\) and \(v\in\ker(P)\), then \(u'v=a'Pv=0\), using symmetry. Thus the two terms are perpendicular, and \(Px\) is the orthogonal projection of \(x\) onto the range. Also, \(x'Px=(Px)'x=(Px)'(Px+(I-P)x)=\|Px\|^2\), as required.

Orthogonality and Rank Additivity

For projection matrices, the algebraic condition \(P_iP_j=0\) exactly expresses orthogonality of their ranges. If the product is zero, then for any vectors \(a,b\), \((P_i a)'(P_j b)=a'P_iP_jb=0\). Conversely, if the ranges are orthogonal, then \(P_iP_jb\) belongs to \(\mathcal{R}(P_i)\), while it is also orthogonal to every vector in \(\mathcal{R}(P_i)\), since \(P_jb\in\mathcal{R}(P_j)\). Hence \(P_iP_jb=0\) for every \(b\), so \(P_iP_j=0\).

Geometric criterion: Suppose \(P_1,\ldots,P_k\) are symmetric idempotent \(n\times n\) matrices with pairwise orthogonal ranges. Then their ranges span \(\mathbb{R}^n\) if and only if \(\sum_{i=1}^k \operatorname{rank}(P_i)=n\). When they span \(\mathbb{R}^n\), every vector has a unique decomposition into one component from each range, and \(\sum_{i=1}^k P_i=I_n\).

Here is a direct proof. Pairwise orthogonal subspaces have a direct sum: if \(\sum_i u_i=0\), with \(u_i\in\mathcal{R}(P_i)\), then taking the squared norm gives \(\|\sum_i u_i\|^2=\sum_i\|u_i\|^2=0\), so every \(u_i=0\). Therefore the dimension of their sum is \(\sum_i\dim\mathcal{R}(P_i)=\sum_i\operatorname{rank}(P_i)\). This sum equals \(n\) precisely when the subspaces span \(\mathbb{R}^n\). In that case each \(x\) has a unique representation \(x=\sum_i u_i\), with \(u_i\in\mathcal{R}(P_i)\). Since the subspaces are orthogonal, \(P_i x=u_i\); summing gives \(\sum_iP_i x=x\) for every \(x\), and hence \(\sum_iP_i=I_n\). The reverse implication follows because a spanning collection of these orthogonal subspaces has dimension \(n\). \(\square\)

Worked Example: Coordinate Subspaces in Three Dimensions

Let \(P_1=\operatorname{diag}(1,1,0)\) and \(P_2=\operatorname{diag}(0,0,1)\). Both matrices are symmetric and idempotent. Their ranges are respectively the first two coordinate directions and the third coordinate direction.

The product is \(P_1P_2=0\), so the ranges are orthogonal. Their ranks are \(2\) and \(1\), and \(2+1=3\), the ambient dimension. Direct addition confirms the identity: \(P_1+P_2=\operatorname{diag}(1,1,1)=I_3\). For \(z=(z_1,z_2,z_3)'\), the two projected squared lengths are \(z'P_1z=z_1^2+z_2^2\) and \(z'P_2z=z_3^2\); their sum is \(z'z\).

Cochran's Theorem in Orthogonal Coordinates

Let \(Z\sim N_n(0,I_n)\). Its distribution is unchanged by an orthogonal change of coordinates. We use this fact to prove the Gaussian conclusion once the projection subspaces form a complete orthogonal decomposition.

Theorem (Cochran's theorem, orthogonal projection form): Let \(Z\sim N_n(0,I_n)\), and let \(P_1,\ldots,P_k\) be symmetric idempotent \(n\times n\) matrices with pairwise orthogonal ranges. If \(\sum_i\operatorname{rank}(P_i)=n\), then \(Z'Z=\sum_{i=1}^k Z'P_iZ\), and the random variables \(Z'P_iZ\) are mutually independent with \(Z'P_iZ\sim\chi^2_{r_i}\), where \(r_i=\operatorname{rank}(P_i)\).

Proof. By the geometric criterion, the ranges form an orthogonal decomposition of \(\mathbb{R}^n\), and \(\sum_iP_i=I_n\). Choose an orthonormal basis for each range. Put the basis vectors together, in order, as the columns of an \(n\times n\) matrix \(U\). Orthogonality between ranges and orthonormality within each basis make \(U\) orthogonal, so \(U'U=I_n\). The columns associated with range \(i\) form a block of \(r_i\) columns.

The coordinates \(W=U'Z\) have the standard \(n\)-variate normal distribution. Indeed, the density of \(Z\) is proportional to \(\exp(-z'z/2)\); the orthogonal transformation preserves \(z'z\) and has Jacobian with absolute value one. Thus the density of \(W\) is \(\prod_{j=1}^n (2\pi)^{-1/2}\exp(-w_j^2/2)\). Its coordinates are independent standard normal random variables.

In these coordinates, \(P_i\) keeps exactly the coordinates in block \(i\) and sets the others to zero. Therefore \(Z'P_iZ\) is the sum of the squares of the \(r_i\) standard normal coordinates in that block. By the definition of the chi-square distribution, \(Z'P_iZ\sim\chi^2_{r_i}\). Distinct blocks use disjoint sets of independent coordinates, so the resulting sums of squares are mutually independent. Finally, summing the projections gives \(\sum_i Z'P_iZ=Z'(\sum_iP_i)Z=Z'Z\). \(\square\)

If \(X\sim N_n(0,\sigma^2I_n)\) with \(\sigma^2>0\), apply the theorem to \(Z=X/\sigma\). Then \(X'P_iX/\sigma^2\sim\chi^2_{r_i}\), the components are mutually independent, and their sum is \(X'X/\sigma^2\).

Worked Example: Two Rotated Coordinate Directions

In \(\mathbb{R}^2\), take the unit vectors \(u=(1,1)'/\sqrt{2}\) and \(v=(1,-1)'/\sqrt{2}\). They are orthogonal because \(u'v=(1-1)/2=0\). The projections onto their spans are \(P_1=uu'=\frac12\begin{pmatrix}1&1\\1&1\end{pmatrix}\) and \(P_2=vv'=\frac12\begin{pmatrix}1&-1\\-1&1\end{pmatrix}\).

Each is symmetric and idempotent: for example, \(P_1^2=u(u'u)u'=uu'=P_1\), and the same calculation applies to \(P_2\). Each has rank \(1\), so the ranks sum to \(2\); also, direct addition gives \(P_1+P_2=I_2\).

For \(Z=(Z_1,Z_2)'\sim N_2(0,I_2)\), \(Z'P_1Z=(Z_1+Z_2)^2/2\) and \(Z'P_2Z=(Z_1-Z_2)^2/2\). The transformed coordinates \((Z_1+Z_2)/\sqrt{2}\) and \((Z_1-Z_2)/\sqrt{2}\) are independent standard normals because the transformation is orthogonal. Thus each projected sum of squares is \(\chi^2_1\), the two are independent, and they sum to \(Z_1^2+Z_2^2\).

When the Matrix Sum Is Given

Sometimes the decomposition is specified algebraically by \(\sum_iP_i=I_n\), rather than by orthogonality and a rank calculation. For symmetric idempotent matrices, that identity itself forces the needed geometry.

Result: If symmetric idempotent matrices \(P_1,\ldots,P_k\) satisfy \(\sum_iP_i=I_n\), then their ranges are pairwise orthogonal and their ranks sum to \(n\).

Proof. Fix an index \(i\) and take any \(x\in\mathcal{R}(P_i)\). Then \(P_ix=x\). From the matrix-sum identity, \[ \|x\|^2=x'x=x'\left(\sum_jP_j\right)x=\sum_j x'P_jx=\sum_j\|P_jx\|^2. \] The last equality uses the projection identity \(x'P_jx=\|P_jx\|^2\). The \(i\)th term is \(\|P_ix\|^2=\|x\|^2\); all the other terms are nonnegative. Consequently every \(j\ne i\) has \(\|P_jx\|^2=0\), so \(P_jx=0\). In particular, \(P_jP_i=0\). Taking transposes gives \(P_iP_j=0\), which means the ranges are orthogonal. Taking traces of the matrix-sum identity now gives \(n=\operatorname{tr}(I_n)=\sum_i\operatorname{tr}(P_i)=\sum_i\operatorname{rank}(P_i)\). The final equality holds because an idempotent matrix has eigenvalues zero or one, so its trace equals its rank. \(\square\)

Worked Example: Checking the Matrix-Sum Condition

Use the two projections from the rotated-direction example. Their sum is \(I_2\), so the result just proved guarantees orthogonal ranges and ranks adding to \(2\). These facts can also be checked directly: \(P_1P_2=uv'v u'=0\), since \(v'u=0\), and \(\operatorname{rank}(P_1)+\operatorname{rank}(P_2)=1+1=2\).

For \(Z\sim N_2(0,I_2)\), the matrix-sum identity gives \(Z'P_1Z+Z'P_2Z=Z'Z\). The coordinate calculation above then establishes that the two terms are independent \(\chi^2_1\) variables. The matrix identity supplies the complete decomposition; Gaussianity supplies the particular distributions and independence.

What Rank Additivity Does—and Does Not—Say

The rank condition is a dimension count, not a substitute for orthogonality. When the ranges are already known to be mutually orthogonal, their dimensions add, and a total of \(n\) dimensions means they exhaust the space. Without orthogonality, ranks can add to \(n\) while the subspaces overlap and leave directions unrepresented.

Worked Example: Equal Ranks but Overlapping Ranges

In \(\mathbb{R}^2\), let \(P_1=P_2=\begin{pmatrix}1&0\\0&0\end{pmatrix}\). Both matrices are symmetric and idempotent, and each has rank \(1\), so their ranks sum to the ambient dimension \(2\). But the ranges are the same line, not orthogonal, and \(P_1+P_2=\begin{pmatrix}2&0\\0&0\end{pmatrix}\ne I_2\).

If \(Z=(Z_1,Z_2)'\sim N_2(0,I_2)\), then \(Z'P_1Z=Z_1^2=Z'P_2Z\). The two quadratic forms are equal, so they cannot be independent: for example, the event \(Z'P_1Z=0\) is exactly the event \(Z'P_2Z=0\), and this event has probability zero but does not establish independence; more directly, their covariance is \(\operatorname{Var}(Z_1^2)=2\), which is not zero. They also fail to sum to \(Z'Z\), since their sum is \(2Z_1^2\), not \(Z_1^2+Z_2^2\).

Common mistake. “The ranks add to \(n\), so the projections must be orthogonal and their quadratic forms independent.” This skips a necessary premise. Rank additivity guarantees exhaustion only for a collection of already orthogonal subspaces. The overlapping-line example has ranks summing to \(n\), but neither the total-sum identity nor independence holds. First verify the projection and orthogonality conditions, or verify the matrix-sum identity; then use the rank count to establish completeness when appropriate.

Key takeaway: Symmetric idempotent matrices are orthogonal projections. Pairwise orthogonal ranges have total dimension \(n\) exactly when they fill \(\mathbb{R}^n\); in that case their projections sum to \(I_n\). Under a centered spherical Gaussian model, orthogonal coordinates turn the projected squared lengths into independent chi-square variables.

Check Your Understanding

Use the projection and Gaussian arguments developed above to answer these questions.

Why does \(x'Px=\|Px\|^2\) hold when \(P\) is symmetric and idempotent?
Suppose three mutually orthogonal projection ranges have dimensions \(2\), \(1\), and \(3\) in \(\mathbb{R}^6\). What follows about their sum and the corresponding Gaussian quadratic forms?
Explain why orthogonality of the ranges implies \(P_iP_j=0\).
Can ranks summing to \(n\) alone establish that a set of projection matrices sums to \(I_n\)? Give the relevant condition that is missing.
Why are the projected sums of squares independent after choosing an orthonormal basis adapted to the ranges?