Common Errors in Matrix-Valued Regression

Why Matrix Bookkeeping Matters

In multivariate regression, a transpose or a divisor that looks harmless can change the meaning of a calculation. A matrix may contain sums of cross-products, estimated covariances, or fitted responses; these objects can have the same dimensions while answering different questions. Careful bookkeeping makes those distinctions explicit.

Use \(n\) for the number of cases, \(p\) for the number of columns in the design matrix, and \(q\) for the number of responses. Write the model as \(Y=XB+U\), where \(Y\) is \(n\times q\), \(X\) is \(n\times p\), \(B\) is \(p\times q\), and \(U\) is \(n\times q\). The rows of \(Y\) and \(U\) correspond to cases; their columns correspond to responses.

The distinction between a data matrix and a coefficient matrix is essential: \(Y\) has one row per case, whereas \(B\) has one row per design column. In particular, the residual SSCP is formed by transposing the residual matrix, not by multiplying it by its transpose in the opposite order. The product \(R'R\) is \(q\times q\); \(RR'\) is \(n\times n\) and summarizes a different set of quantities.

Worked Examples: Dimensions, SSCP, and Covariance

Residual Degrees of Freedom and the Covariance Estimate

Let \(r=\operatorname{rank}(X)\), and let \(P_X\) be the orthogonal projection onto the column space of \(X\). Least-squares residuals satisfy \(R=(I-P_X)Y\). The residual projection \(M=I-P_X\) is symmetric and idempotent, and its rank and trace are both \(\nu=n-r\). This remains the appropriate degree-of-freedom count if the design columns are linearly dependent; replacing \(r\) by \(p\) in that situation is an error.

Proof. Since the mean \(XB\) lies in the column space of \(X\), \(M X B=0\), so \(R=MU\). Therefore \(R'R=U'M'U=U'MU\). For response coordinates \(a\) and \(b\), the corresponding entry is \[ (U'MU)_{ab}=\sum_{i=1}^n\sum_{j=1}^n m_{ij}u_{ia}u_{jb}. \] Independence of different rows and their zero means imply \(E[u_{ia}u_{jb}]=0\) when \(i\ne j\), while \(E[u_{ia}u_{ib}]=\Sigma_{ab}\). Taking expectations yields \[ E[(U'MU)_{ab}]=\sum_{i=1}^n m_{ii}\Sigma_{ab} =\operatorname{tr}(M)\Sigma_{ab} =\nu\Sigma_{ab}. \] This holds for each pair \(a,b\), so \(E[R'R]=\nu\Sigma\). Dividing by positive \(\nu\) proves unbiasedness. \(\square\)

The assumptions matter. The formula does not automatically apply when case errors are correlated or have different covariance matrices. Nor does unbiasedness make \(E/\nu\) a universally correct divisor for every matrix summary: it applies to the residual covariance estimate under the stated model assumptions.

A Gaussian Maximum-Likelihood Caveat

A related pitfall is to state without qualification that the Gaussian maximum-likelihood estimate of \(\Sigma\) is \(E/n\). That expression is valid only when the covariance maximum exists in the positive-definite parameter space. In the Gaussian model with the mean fitted by least squares, the covariance-dependent part of the likelihood is proportional to \[ |\Sigma|^{-n/2}\exp\{-\tfrac12\operatorname{tr}(\Sigma^{-1}E)\}. \]

If \(E\) is positive definite, maximizing this likelihood over positive-definite \(\Sigma\) gives \(\hat{\Sigma}_{\mathrm{ML}}=E/n\). If \(E\) is singular, there is no finite maximum over positive-definite covariance matrices. To see why, diagonalize \(E\) in an orthonormal basis, with \(s<q\) positive eigenvalues and at least one zero eigenvalue. Choose \(\Sigma\) diagonal in the same basis, using eigenvalue \(\lambda_j/n\) for each positive eigenvalue \(\lambda_j\) of \(E\), and eigenvalue \(\epsilon>0\) in each zero-eigenvalue direction. The trace term stays fixed as \(\epsilon\) decreases, while \(|\Sigma|^{-n/2}\) grows without bound. Thus the likelihood is unbounded as \(\epsilon\downarrow0\), and no finite positive-definite covariance MLE exists. Since \(\operatorname{rank}(E)\leq\nu\), \(\nu<q\) guarantees singularity, though \(\nu\geq q\) alone does not guarantee that \(E\) is positive definite.

Reading Off-Diagonal Entries and Response Correlation

For residual columns \(r_a\) and \(r_b\), \(E_{ab}=r_a'r_b\). A positive value means the residuals tend to have products with positive sum; a negative value means their products tend to have negative sum. The sign by itself is not a correlation, because the cross-product depends on the scale of both residual columns.

When the diagonal entries are positive, the residual correlation calculated from \(E\) is \[ \frac{E_{ab}}{\sqrt{E_{aa}E_{bb}}}. \] Dividing every entry by the same positive number, such as \(\nu\), leaves this correlation unchanged. But residual response correlation is not generally the same as raw response correlation: raw responses may be associated because their fitted means vary together with the predictors. Residual correlation describes association remaining after the specified mean structure has been fitted.

Common Errors and Their Corrections

• Transposing the wrong product. \(R'R\) is \(q\times q\) and contains cross-products among responses. \(RR'\) is \(n\times n\), with entries pairing cases. Check dimensions and intended interpretation before multiplying.
• Calling the SSCP a covariance matrix. \(E=R'R\) contains sums. Under the stated common-covariance assumptions, \(E/\nu\), not \(E\), is the unbiased residual covariance estimate.
• Using \(n-p\) when \(X\) is rank deficient. Residual degrees of freedom are \(n-r\), with \(r=\operatorname{rank}(X)\). The number of listed columns is not necessarily the number of independent fitted directions.
• Calling \(E/n\) the Gaussian MLE in every case. It is the MLE only when \(E\) is positive definite. If \(E\) is singular, the likelihood over positive-definite covariance matrices is unbounded.
• Confusing residual and raw response correlation. Off-diagonal residual cross-products describe association left after fitting the mean model; they do not automatically describe the marginal association of the observed responses.