Deriving Multivariate Least Squares by Projection

Why Least Squares Is a Projection Problem

In a multivariate linear model, each case has several responses, while the same design matrix describes the predictors for every response. The least-squares problem is therefore to choose one coefficient matrix that fits all response columns jointly, minimizing their total squared error. The trace formulation makes that total precise; the geometry of projection explains why the minimizing fitted matrix has the form it does.

Let $Y$ be an $n\times q$ response matrix and $X$ an $n\times p$ design matrix. The $i$th row of $Y$ contains the $q$ responses for case $i$, and $X$ has one row per case. The coefficient matrix $B$ is $p\times q$, so $XB$ is $n\times q$, matching $Y$. Assume throughout that $X$ has full column rank $p$; equivalently, $X'X$ is invertible. This condition gives the unique coefficient matrix in the formula derived below.

Definition: For a $p\times q$ coefficient matrix $B$, define the least-squares objective by $$ Q(B)=\operatorname{tr}\big((Y-XB)'(Y-XB)\big). $$ For any matrix $A$, its squared Frobenius norm is $\|A\|_F^2=\operatorname{tr}(A'A)$. Thus $Q(B)=\|Y-XB\|_F^2$: it is the sum of the squared entries of the residual matrix.

Writing the objective as a trace does not weight one response column differently from another. If $y_j$ and $b_j$ are column $j$ of $Y$ and $B$, respectively, then $Q(B)=\sum_{j=1}^q \|y_j-Xb_j\|^2$. The coefficient columns are estimated together in matrix notation, but the criterion adds the ordinary squared errors across all cases and responses.

From the Objective to the Coefficient Matrix

Expanding the trace exposes the quadratic function of $B$. The trace identities $\operatorname{tr}(A)=\operatorname{tr}(A')$ and $\operatorname{tr}(CD)=\operatorname{tr}(DC)$, when the products are defined, give

\begin{aligned} Q(B) &=\operatorname{tr}(Y'Y)-\operatorname{tr}(Y'XB) -\operatorname{tr}(B'X'Y)+\operatorname{tr}(B'X'XB)\\ &=\operatorname{tr}(Y'Y)-2\operatorname{tr}(B'X'Y) +\operatorname{tr}(B'X'XB). \end{aligned}

The two middle terms are equal scalars, since one is the transpose of the other. The next result identifies the minimizer and proves it is unique.

Theorem: If $X$ has full column rank, the unique minimizer of $Q(B)$ over all $p\times q$ matrices is $$ \widehat B=(X'X)^{-1}X'Y. $$ It satisfies the normal equations $X'X\widehat B=X'Y$.

Proof. Put $A=X'X$ and $C=X'Y$. Full column rank makes $A$ positive definite: for any nonzero $v\in\mathbb{R}^p$, $v'Av=\|Xv\|^2>0$. Define $B_0=A^{-1}C$. Then $AB_0=C$. For any $B$, let $D=B-B_0$. Substituting $B=B_0+D$ into the expanded objective gives

\begin{aligned} Q(B)-Q(B_0) &=-2\operatorname{tr}(D'C) +\operatorname{tr}\big((B_0+D)'A(B_0+D)-B_0'AB_0\big)\\ &=-2\operatorname{tr}(D'C) +\operatorname{tr}(D'AB_0)+\operatorname{tr}(B_0'AD) +\operatorname{tr}(D'AD). \end{aligned}

Because $AB_0=C$ and $A$ is symmetric, the first two linear terms in the last line cancel: $\operatorname{tr}(D'AB_0)=\operatorname{tr}(D'C)$, and $\operatorname{tr}(B_0'AD)=\operatorname{tr}(C'D)=\operatorname{tr}(D'C)$. Therefore

Q(B)-Q(B_0)=\operatorname{tr}(D'AD).

Writing $D$ by columns $d_1,\ldots,d_q$ shows that $\operatorname{tr}(D'AD)=\sum_{j=1}^q d_j'Ad_j\geq 0$, with equality only if every $d_j=0$, since $A$ is positive definite. Thus $B_0$ is the unique minimizer. Its definition gives $B_0=(X'X)^{-1}X'Y$, and multiplying by $X'X$ gives the normal equations. $\square$

Worked Example: Two Responses With an Intercept and a Slope

Suppose three cases have predictor values $0,1,2$, and two measured responses are recorded in

X=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}, \qquad Y=\begin{pmatrix}1&0\\2&1\\2&4\end{pmatrix}.

The first column of $X$ is the intercept. Calculate

X'X=\begin{pmatrix}3&3\\3&5\end{pmatrix}, \qquad (X'X)^{-1}=\frac{1}{6}\begin{pmatrix}5&-3\\-3&3\end{pmatrix}, \qquad X'Y=\begin{pmatrix}5&5\\6&9\end{pmatrix}.

For example, the inverse is correct because multiplying the two matrices $X'X$ and $(X'X)^{-1}$ gives the $2\times2$ identity matrix. Thus

\widehat B=(X'X)^{-1}X'Y =\frac{1}{6}\begin{pmatrix}5&-3\\-3&3\end{pmatrix} \begin{pmatrix}5&5\\6&9\end{pmatrix} =\begin{pmatrix}7/6&-1/3\\1/2&2\end{pmatrix}.

The first column gives the intercept and slope for response 1; the second gives them for response 2. The fitted response matrix is

X\widehat B= \begin{pmatrix}7/6&-1/3\\5/3&5/3\\13/6&11/3\end{pmatrix}.

The residuals are $Y-X\widehat B$, with columns $(-1/6,1/3,-1/6)'$ and $(1/3,-2/3,1/3)'$. Their squared sums are $1/6$ and $2/3$, so the minimized objective is $Q(\widehat B)=1/6+2/3=5/6$. This calculation illustrates that both responses use the same design matrix, while each has its own coefficient column.

Why the Formula Is a Projection

The columns of $X$ span a subspace of $\mathbb{R}^n$, called the design space. Every column of $XB$ lies in this space. The matrix

P_X=X(X'X)^{-1}X'

maps each vector in $\mathbb{R}^n$ to its orthogonal projection onto the design space. To check this, $P_X'=P_X$, because $X'X$ and its inverse are symmetric. Also, $P_X^2=X(X'X)^{-1}X'X(X'X)^{-1}X'=P_X$. Its columns lie in the column space of $X$, and $P_X X=X$, so it leaves every vector in that space unchanged. These properties characterize the orthogonal projection onto that space.

Theorem: Under the full-column-rank condition, the least-squares fitted matrix is the projection of $Y$ onto the design space, column by column: $$ X\widehat B=P_XY. $$ Among all matrices $XB$ with columns in the design space, this is the one minimizing the total squared distance to $Y$.

Proof. The coefficient formula immediately gives $X\widehat B=X(X'X)^{-1}X'Y=P_XY$. To verify the minimizing property, take any $B$ and split the difference into two matrices:

$$ Y-XB=(I-P_X)Y+(P_XY-XB). $$

Every column of $(I-P_X)Y$ is orthogonal to the design space: $X'(I-P_X)=X'-X'X(X'X)^{-1}X'=0$. Every column of $P_XY-XB$ lies in the design space, because both $P_XY$ and $XB$ do. The two terms are therefore orthogonal column by column. Expanding the squared Frobenius norm, their cross-products are zero, giving

Q(B)=\|(I-P_X)Y\|_F^2+\|P_XY-XB\|_F^2.

The first term does not depend on $B$, and the second is nonnegative. It equals zero when $XB=P_XY$, which is attained at $B=\widehat B$. Full column rank ensures that no other $B$ gives the same fitted matrix: if $X(B-\widehat B)=0$, each column of $B-\widehat B$ must be zero. Hence the minimizer is unique. $\square$

Worked Example: A Group-Indicator Design

Consider four cases split into two groups, with a separate group mean for each of two responses. The design and response matrices are

X=\begin{pmatrix}1&0\\1&0\\0&1\\0&1\end{pmatrix}, \qquad Y=\begin{pmatrix}2&0\\4&2\\6&1\\8&3\end{pmatrix}.

Here $X'X=\begin{pmatrix}2&0\\0&2\end{pmatrix}$, so $(X'X)^{-1}=\frac12 I$. Also,

X'Y=\begin{pmatrix}6&2\\14&4\end{pmatrix}, \qquad \widehat B=\frac12 X'Y =\begin{pmatrix}3&1\\7&2\end{pmatrix}.

The fitted rows are $(3,1),(3,1),(7,2),(7,2)$: each row is the mean response vector for its group. The residual rows are $(-1,-1),(1,1),(-1,-1),(1,1)$. Each row has squared length $2$, so $Q(\widehat B)=4(2)=8$. The projection interpretation is visible directly: each response column is replaced within each group by its group average, and the deviations from that average sum to zero within the group.

Worked Example: Checking the Completed-Square Increase

Let $X=(1,0,1)'$ and let $Y=( (1,2)',(3,0)',(5,4)' )'$, so $Y$ has rows $(1,2),(3,0),(5,4)$. Then $X'X=2$, $X'Y=(6,6)$, and $\widehat B=(3,3)$. The fitted rows are $(3,3),(0,0),(3,3)$. The residual rows are $(-2,-1),(3,0),(2,1)$, giving $Q(\widehat B)=5+9+5=19$.

Now use the alternative coefficient row $B=(2,4)$. Its fitted rows are $(2,4),(0,0),(2,4)$, and its residual rows are $(-1,-2),(3,0),(3,0)$. Therefore $Q(B)=5+9+9=23$. The completed-square identity predicts an increase of $(B-\widehat B)'(X'X)(B-\widehat B)$. Since $B-\widehat B=(-1,1)$, that increase is $2[(-1)^2+1^2]=4$, agreeing with $23-19=4$.

Interpretation and Common Pitfalls

Projection explains the central role of the design matrix: least squares can choose only fitted columns within the design space. It also explains why the same $X$ yields the same projection operation for every response, even though the estimated coefficients differ across responses. The algebraic formula is useful for direct calculation; the projection view clarifies what the calculation is doing.

Using an inverse of $X$. In general, $X$ is rectangular, so $X^{-1}$ is not defined. The formula uses the inverse of the $p\times p$ matrix $X'X$, and that inverse requires full column rank.
Forgetting the responses are columns. $B$ has $q$ columns, not one coefficient vector for the entire matrix $Y$. Each response has its own coefficient vector, arranged as a column of $B$.
Assuming the trace minimizes a different criterion. Here the trace is exactly the sum of squared residual entries. It does not introduce a covariance weighting between responses.
Thinking any coefficient matrix can fit $Y$ exactly. The fitted columns must lie in the design space. If a response column is outside that space, its projection leaves a nonzero discrepancy.

A common incorrect argument is, “Because the responses are correlated, least squares must change the coefficient estimate by weighting the columns according to their covariance.” That does not follow from this objective: $Q(B)$ is the unweighted sum of squared entries. Under the stated full-rank condition, its minimizer is the same matrix $(X'X)^{-1}X'Y$, regardless of whether a covariance model for the responses is specified.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

What are the dimensions of $X$, $Y$, and $B$, and why must $XB$ have the same dimensions as $Y$?
State the full-rank condition needed for the inverse $(X'X)^{-1}$ in the least-squares formula.
Explain why $\operatorname{tr}((Y-XB)'(Y-XB))$ is the sum of squared residual entries.
Why does the completed-square identity show that the stated coefficient matrix is the unique minimizer?
In the projection argument, why are $(I-P_X)Y$ and $P_XY-XB$ orthogonal?

Deriving Multivariate Least Squares by Projection

What You'll Learn